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A mathematical  model  for  the  dynamic  analysis  of  the  riser  system  with  the  inclu- 
sion of  internal  flow  and  nonlinear  effects  due  to  large  structural  displacements  is  devel- 
oped to  investigate  the  effect  of  internal  flow  on  marine  riser  dynamics.  The  riser  system 
accounts  for  the  nonlinear  boundary  conditions  and  includes  a steady  flow  inside  the  pipe 
which  is  modeled  as  an  extensible  or  inextensible,  tubular  beam  subject  to  nonlinear 
three-dimensional  hydrodynamic  loads  such  as  current  or  wave  excitation.  Combining 
Iwan-Blevin's  model  into  the  approximated  form  of  the  nonlinear  model  derived  in  this 
study,  current-vortex  model  is  derived  to  examine  the  effect  of  internal  flow  on  vortex- 
induced  vibration  due  to  inline  current. 

A semi-analytical  method  using  series  expansion  is  employed  to  derive  Eigenvalue 
problem  to  facilitate  the  evaluation  of  the  system  frequencies,  and  its  validity  is  given 
through  the  comparison  of  the  solutions  with  the  conventional  method  using  system  ma- 
trices. Galerkin’s  finite  element  approximation  and  time-incremental  operator  are  im- 
plemented to  derive  the  matrix  equation  of  equilibrium  for  the  finite  element  system  and 
the  extensibility  or  inextensibility  condition  is  used  to  reduce  the  degrees  of  freedom  of 
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the  system  and  the  required  computational  time  in  the  case  of  a nonlinear  model.  The 
algorithm  is  implemented  to  develop  computer  programs  used  in  several  numerical 
applications. 

The  investigations  of  the  effect  of  internal  flow  on  system  frequency,  vortex- 
induced  vibration,  and  vibration  due  to  current  or  wave  loading  are  performed  according 
to  the  change  of  various  parameters  such  as  top  tension,  internal  flow  velocity,  current 
velocity,  wave  period,  and  so  on.  It  is  found  that  the  effect  of  internal  flow  can  be  con- 
trolled by  the  increase  of  top  tension.  However,  careful  consideration  has  to  be  given  in 
the  design  point,  particularly  for  the  long  riser  under  the  harmonic  loading  such  as 
vortex- shedding  or  waves.  And  also,  the  consideration  of  nonlinear  effects  due  to  large 
structural  displacements  increases  the  effect  of  internal  flow  on  riser  dynamics. 
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CHAPTER  1 
INTRODUCTION 


1.1  General  Description 

The  marine  riser  is  a conductor  pipe  used  in  floating  drilling  operations  to  convey 
drilling  fluid  and  to  guide  tools  between  the  drilling  vessel  and  the  well  head  at  the  ocean 
floor.  The  riser  can  be  a single  rigid  pipeline  for  the  drilling  riser  system,  or  a bundle  of 
flow  lines  that  are  assembled  as  an  integral  unit  for  the  production  riser  system.  At  the 
surface,  the  riser  system  is  connected  to  a surface  vessel  such  as  the  tension  leg  platform 
(TLP),  floating  production  and  storage  unit  (FPSU).  At  the  bottom,  it  may  be  hanging 
free  or  connected  to  a blow-out  prevention  (BOP)  stack. 

Risers  are  often  used  in  the  conventional  sense  to  resist  the  applied  forces  and  to 
transmit  loads  in  the  surface  structure  by  its  bending,  shear  and  axial  stiffness.  To  avoid 
the  collapse  and  buckle  of  the  riser  by  its  own  weight  and  to  reduce  the  excess  load  on 
the  manifold  in  the  deep  sea,  a constant  top  tension  is  often  applied.  Buoyant  material  is 
strapped  around  the  riser  at  intervals  to  provide  an  additional  buoyancy  and  to  reduce  the 
required  top  tension.  Various  end-conditions  for  a riser  system  are  possible,  for  example, 
hinged-hinged,  fixed-free,  etc.  These  are  some  of  the  features  of  the  marine  riser  system. 

The  offshore  industry  is  moving  into  deeper  waters  and  more  hostile  environments. 
When  exploratory  drilling  operations  move  into  deeper  water,  the  dynamic  behaviors  of 
riser  should  be  more  carefully  considered  because  of  the  long  free  span.  There  are  many 
factors  that  influence  on  the  behavior  of  a riser,  such  as  the  offshore  environment  (wave 
and  current),  the  behaviors  of  platform  (movements  or  mean  offsets),  and  the  properties 
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of  the  riser  itself  (dimensions,  weight,  top  tension,  and  internal  flow).  Most  of  these  fac- 
tors have  been  considered  in  earlier  studies  and  many  related  papers  have  been  published. 
However,  it  was  not  until  the  publication  of  a paper  (Moe  and  Chucheepsakul,  1988)  that 
the  effect  of  internal  flow  on  riser  dynamics  was  investigated.  Before  Moe's  investigation 
it  was  already  known  that  the  internal  flow  might  cause  a dynamic  instability  or  buckling 
of  a horizontal  pipeline  supported  above  ground  and  such  effect  has  been  considered  in 
design  and  analysis  of  pipeline.  Since  the  dynamics  of  pipe  conveying  fluid  has  been  stu- 
died fairly  extensively  over  the  past  40  years,  the  basic  techniques  were  recently  adopted 
for  the  dynamic  analysis  of  a marine  riser  with  internal  flowing  fluid. 

Figure  1.1  is  a schematic  of  a marine  riser  conductor  that  contains  internal  fluid 
flow  and  is  subject  to  environmental  forces.  When  the  internal  fluid  travels  inside  the 
curved  path  along  the  deflected  riser,  it  experiences  centrifugal  and  coriolis  accelerations 
due,  respectively,  to  the  curvature  of  the  riser  and  the  relative  motion  of  fluid  to  time  de- 
pendent riser  motion.  Those  accelerations  exert  against  the  riser  which,  in  turn,  affect  the 
dynamic  behavior  of  the  riser  and  cause  riser  vibrations.  In  addition,  the  riser  vibration 
could  be  "locked"  with  the  flow-induced  vortex  shedding  such  that  they  resonate  together 
to  produce  large  deflections  and  stresses.  This  phenomenon  could  lead  to  the  failure  of 
the  riser  system  prematurely. 

In  addition  to  this  interaction  problem,  the  riser  behavior  is  also  nonlinear.  The 
nonlinearities  are  mainly  of  two  origins  that  result  from  flow-induced  drag  force  and 
from  geometric  nonlinearities  due  to  large  structural  deflections.  The  latter  becomes  par- 
ticularly significant  for  long  risers  over  1000  m.  Such  nonlinearities  are  due  to  large 
deflections  and  slopes,  three-dimensional  bending,  extensibility,  torsion,  dependency  of 
the  hydrodynamic  loads  on  the  riser  deformation  and  nonlinear  boundary  conditions. 
Three-dimensional  bending  is  also  important  for  such  risers  because  the  hydrodynamic 
loading  due  to  the  surface  waves,  the  ocean  current,  the  relative  motion  of  the  riser  with 
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Figure  1.1:  Configuration  of  Riser  System  with  the  Inclusionn  of  Internal  Flow 

respect  to  water  particle  and  the  position  of  the  supporting  offshore  platform,  is  actually 
three-dimensional.  Extensibility  effects  are  important  for  production  risers  made  of  flex- 
ible composite  material.  Torsional  couples,  concentrated  at  the  upper  end  of  the  riser, 
arise  from  the  variation  of  the  relative  orientation  of  the  supporting  offshore  platform 
and  the  riser.  The  hydrodynamic  loading  exerted  on  the  riser  depends  on  the  orientation 
of  the  riser  tubes  with  respect  to  the  surrounding  flow  and  is  deformation  dependent.  The 
displacements  resulting  from  these  nonlinear  effects,  in  turn,  modify  the  particle  motion 
of  the  internal  fluid.  The  system  becomes  a complicated  interaction  problem. 
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In  summary,  the  primary  objectives  of  this  study  are  1)  to  develop  a mathematical 
model  for  the  analysis  of  the  riser  system  with  the  inclusion  of  internal  flow  and  the 
aforementioned  nonlinear  effects,  2)  to  solve  the  model  numerically,  and  3)  to  examine 
the  effect  of  the  internal  flow  on  riser  dynamics.  In  this  study,  Iwan-Blevin's  model 
(1990),  which  is  one  of  the  self-excited  oscillator  models,  will  be  implemented  as  a vor- 
tex shedding  model  to  examine  the  effect  of  the  internal  flow.  Semi-analytical  solution 
using  series  expansions  and  numerical  methods  using  system  matrices  are  both  employed 
to  derive  the  Eigenvalue  problem  to  facilitate  the  evaluation  of  system's  natural  frequen- 
cies. For  the  solution  of  forced  vibration,  Galerkin's  finite  element  approximation  is  im- 
plemented and  the  Blevins-Iwan  model  (1990)  is  adapted  to  solve  vortex-induced 
vibration  due  to  inline  current.  The  effect  of  internal  flow  on  system  frequencies,  vortex- 
induced  vibration,  wave-induced  vibration,  internal  stresses  as  well  as  the  kinematics  of 
riser  are  investigated  according  to  the  change  of  internal  flow  or  current  velocity,  water 
depth,  top  tension,  wave  period  and  so  on. 

1.2  Literature  Review 

The  literature  review  contains  three  categories:  1)  Instability  of  fluid-conveying 
pipe  in  which  the  basic  concept  is  presented,  2)  The  effect  of  internal  flow  on  riser  dy- 
namics in  which  the  previous  researches  of  this  subject  are  discussed,  and  3)  Riser  analy- 
sis in  which  rational  ideas  for  extending  previous  works  are  given. 

1.2.1  Instability  of  Fluid-Conveving  Pipe 

The  vibration  of  pipe  conveying  fluid  did  not  attract  much  attention  until  1950 
when  severe  vibrations  of  the  Trans-Arabian  pipeline  were  observed.  Ashley  and  Havi- 
land  (1950)  studied  the  bending  vibrations  of  a pipeline  in  conjunction  with  the  flow- 
induced  vibrations  of  the  Trans-Arabian  pipeline.  Housner  (1952)  was  the  first  to 
correctly  derive  the  governing  equation  of  motion  and  predict  instability.  He  concluded 
that  the  internal  flow  velocity  might  cause  a dynamic  coupling  of  the  simple  modes  of 
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vibration  so  that  the  normal  modes  of  vibration  were  of  complex  shape  with  90-  degree 
out  of  phase  components.  At  low  fluid  velocities  there  was  negligible  effect  on  the  vibra- 
tion of  the  pipeline,  but  at  a certain  high  critical  velocity  the  fluid  flow  caused  a dynamic 
instability.  The  induced  centrifugal  and  coriolis  accelerations  become  important  in  the 
deflected  pipeline. 

The  type  of  instability  depends  on  the  end  conditions  of  the  pipe.  Pipes  supported 
at  both  ends  blow  out  and  buckle  when  the  flow  velocity  exceeds  the  critical  velocity 
(Housner,  1952;  Homes,  1978).  Benjamin  (1961),  dealing  with  the  general  dynamic 
problem  of  articulated  pipe  systems  conveying  fluid,  on  the  other  hand,  showed  that 
when  such  systems  possess  one  free  end,  unstable  oscillations  may  develop  for  sufficient- 
ly high  velocities.  Paidoussis  (1966)  demonstrated  the  existence  of  unstable  oscillations 
of  continuously  flexible  tubular  cantilevers  conveying  fluid.  In  his  study,  a cylindrical 
body  was  immersed  in  fluid  flowing  parallel  to  the  direction  of  internal  flow  and  the 
forces  exerted  by  the  fluid,  except  those  due  to  fluid  friction,  were  caused  from  lateral 
acceleration.  Later,  he  (Paidoussis,  1970)  showed  that  cantilever  pipes  on  the  ground 
were  flailed  with  large  amplitude  once  the  internal  flow  velocity  exceeds  a certain  critical 
velocity. 

Chen  (1973)  presented  the  out-of-plane  motion  of  uniformly  curved  tubes  contain- 
ing flowing  fluid.  Hamilton's  principle  was  used  to  derive  both  in-plane  and  out-of-plane 
equations  of  motion  and  a general  solution  for  the  natural  frequency  was  obtained  ana- 
lytically. Hill  and  Davis  (1974)  dealt  with  the  effect  of  initial  forces,  which  could  arise  in 
a clamped  tube  as  it  is  loaded  by  a pressurized  flowing  fluid  on  an  initially  planar  curved 
tube.  They  obtained  numerical  solutions  by  using  the  finite  element  technique. 

Paidoussis  and  Issid  (1974)  investigated  the  dynamics  and  instability  of  flexible 
pipes  containing  flowing  fluid,  where  the  flow  velocity  is  either  entirely  constant,  or  with 
a small  harmonic  component  superposed.  Paidoussis  (1986)  re-examined  the  dynamics 
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and  stability  of  short  tubes  by  means  of  Timoshenko  beam  theory  and  a three- 
dimensional  fluid  model  for  the  flow,  rather  than  the  plug-flow  model  commonly  uti- 
lized. In  their  study,  it  was  shown  that  for  long  tubes,  Euler-Bemoulli  beam  theory  and  a 
plug-flow  model  were  perfectly  adequate. 

In  some  of  the  latest  works,  Paidoussis  (1988)  and  Blevins  (1990)  provided  general 
reviews  on  dynamics  of  fluid-conveying  pipes  which  listed  many  references  about  this 
subject. 

1.2.2  The  Effect  of  Internal  Flow  on  Riser  Dynamics 

Although  several  studies  have  investigated  the  vibration  of  a pipeline  conveying 
fluid  supported  above  ground,  there  are  only  a few  papers  dealing  with  the  effect  of  inter- 
nal flow  on  marine  riser  dynamics. 

Moe  (1988)  was  the  first  who  considered  the  forces  due  to  the  internal  flowing  flu- 
id as  a dynamically  forcing  component  acting  on  the  interior  wall  of  the  riser  and  derived 
the  governing  equation  of  motion.  Before  that,  the  loading  due  to  the  internal  fluid  was 
included  in  the  internal  tension  of  the  riser  as  a fluid  static  force,  that  is,  centrifugal  and 
coriolis  acceleration  due  to  the  internal  flow  were  largely  negligible. 

Wu  and  Lou  (1991)  developed  a mathematical  model  for  the  lateral  bending  vibra- 
tion of  a marine  riser  and  examined  the  effect  of  internal  flow  and  bending  rigidity  of  the 
pipe  on  the  dynamic  behavior  of  the  riser.  Their  mathematical  model  included  the  steady 
flow  inside  the  pipe  together  with  other  factors  such  as  currents  or  wave  excitation  and 
used  a singular  perturbation  technique  to  solve  it  under  the  condition  that  the  rigidity  pa- 
rameter, £ = EI/L3 , is  sufficiently  small  for  deep  water  risers.  It  was  found  from  their 
results  that  the  internal  flow  reduced  the  effect  of  top  tension,  but  the  riser  motion  was 
not  significantly  affected  when  the  top  tension  of  the  riser  was  relatively  high.  However, 
the  problem  deserves  further  investigation  since  the  system  solved  was  linear  and  the  per- 
turbation technique  was  not  valid  for  low  top  tension  cases. 
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Chen  (1992)  derived  the  governing  equation  with  the  inclusion  of  dynamic  force 
for  lateral  vibration  by  applying  Hamilton's  principle.  The  natural  frequencies  and  the 
mode  shapes  were  formulated  and  presented.  The  critical  buckling  and  significant  veloci- 
ties which  associate  the  internal  flowing  fluid  with  system  integrity  were  also  presented. 
In  his  paper,  it  was  concluded  that  the  natural  frequencies  could  be  reduced  drastically 
by  the  fluid  dynamic  force  if  the  tension  was  insufficient  for  a neutrally  buoyant  riser 
system,  and  that  the  dynamic  force  had  less  influence  on  a positively  or  negatively  buoy- 
ant riser  system.  However,  his  system  is  also  a linear  model  as  was  that  of  Wu  and  the 
method  employed  by  him  to  estimate  the  natural  frequencies  and  the  modal  shapes  ap- 
peared to  be  quite  approximate. 

1.2.3  Riser  Analysis 

For  the  riser  analysis  without  considering  the  effect  of  internal  flow,  considerable 
efforts  have  been  made  over  the  years.  A marine  riser  can  be  treated  as  a long,  slender 
tubular  beam-column  (Bemoulli-Euler  type)  connected  at  upper  and  lower  ball  joint.  For 
linear  analysis,  several  methods,  e.g.,  finite  difference,  finite  element  method,  or  analyti- 
cal method,  etc.,  were  used.  Burke  (1974)  studied  the  dynamic  riser  response  and  dis- 
cussed the  effects  of  increasing  water  depth  on  riser  design  and  operations.  Dareing  and 
Huang  (1979)  explored  the  application  of  modal  analysis  for  calculating  marine  riser 
time-dependent  stresses  and  suggested  the  simplification  of  dynamic  response  calcula- 
tions, using  only  a limited  number  of  eigenvalues  and  eigenfunctions.  Chakrabati  and 
Frampton  (1982)  presented  a review  on  riser  analysis  techniques. 

With  increasing  water  depth,  the  nonlinear  behaviors  of  the  riser  become  more 
prominent  as  mentioned  in  section  1.1.  The  study  of  riser  analysis  including  nonlineari- 
ties has  been  severely  limited  by  the  complex  governing  equations  describing  its  equilib- 
rium configuration.  With  the  advances  of  numerical  methods,  however,  it  became 
feasible  to  solve  nonlinear  systems.  The  equations  for  large  amplitude  three-dimensional 
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inextensional  motions  of  beams  were  derived  by  Nordgren  (1974),  with  the  assumption 
of  constant  principal  moments  of  inertia,  negligible  rotatory  inertia,  and  the  uncoupled 
torsional  motion.  Bruce  and  Michael  (1977)  described  a mathematical  model  and  solu- 
tion technique  for  a system  with  coupled  dynamic  axial  and  lateral  responses  of  a riser 
column. 

Fellipa  and  Chung  (1981)  formulated  a mathematical  model  and  implemented  a 
finite  element  method  for  solving  nonlinear  static  equilibrium  configurations  of  deep  wa- 
ter risers.  Also,  a transient  analysis  was  developed  by  Chung  (1981)  for  the  determina- 
tion of  nonlinear  motion  by  considering  nonlinear  static  configuration  as  an  initial 
condition  for  dynamic  analysis.  A tribute  to  their  work  was  the  consistent  treatment  of 
hydrodynamic  loading. 

Garrett  (1982)  presented  a three-dimensional  finite  element  model  of  an  inextens- 
ible  elastic  rod  with  equal  principal  stiffness.  The  model  permitted  large  deflections  and 
finite  rotations  and  accounted  for  tension  variation  along  its  length.  Its  use  in  static  analy- 
sis was  described  and  time  integration  method  for  dynamic  analysis  was  developed. 

Konuk  (1982)  provided  a general  foundation  for  developing  rigorous  formulation 
of  problems  involving  marine  pipelines  with  twist.  A continuation  technique  was  pro- 
posed to  handle  convergence  problems. 

Safai  (1983)  developed  a method  for  automatically  updating  the  structural  geome- 
try during  the  dynamic  analysis  for  a system  in  which  the  bending,  axial  and  torsional 
vibrations  are  uncoupled.  The  advantage  of  his  method  was  that  it  could  be  applied  to 
risers  in  great  water  depth  with  relatively  large  top  displacements,  for  which  the  global 
rotation  of  each  riser  element  easily  exceeded  the  limit  considered  to  discriminate  large 
and  small  deformations,  i.e.  8 degree. 

Kim  and  Triantafyllou  (1984)  studied  the  nonlinear  dynamics  of  long,  slender  cyl- 
inders assuming  moderately  large  deformations  and  no  longitudinal  excitation  modes.  In 
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their  paper,  the  nonlinearity  in  axial  strain  was  represented  through  variable  tension  de- 
pending on  time  and  space,  while  their  bending  equations  were  linear  and  uncoupled. 

McNamara  and  Lane  (1984)  presented  an  efficient  method  based  on  the  finite  ele- 
ment approach  using  convected  coordinates  for  arbitrary  large  rotations  and  established 
an  associated  computer  program  for  computing  the  gross  motions  and  internal  forces  of 
offshore  columns  and  loading  arms. 

Huang  and  Chucheepsakul  (1985)  introduced  a method  of  static  analysis  for  risers 
experiencing  large  displacements.  The  method  represented  a modified  functional  involv- 
ing multipliers  to  account  for  the  arc  length  being  unvaried  with  neglected  tension,  and 
used  the  exact  curvature  of  radius  in  their  formulation.  This  method  is  suitable  for  ana- 
lyzing a riser  having  a known  top  tension  and  a possible  slippage  at  the  top  slip  joint. 

Kokarakis  and  Bemitsas  (1985)  formulated  the  problem  of  static  three-dimensional 
nonlinear,  large  deformation  response  of  a riser  within  small  strain  theory  and  then 
solved  it  numerically  by  using  an  incremental  finite  element  algorithm,  which  involved  a 
prediction-correction  scheme.  The  extensibility  or  inextensibility  condition  was  used  as 
the  constitutive  relation  in  the  tangential  direction.  Torsion  and  bending  were  coupled. 
The  external  loads  and  boundary  conditions  were  deformation  dependent.  Later,  they 
(Bemitsas  and  Kokarakis,  1987)  employed  a previously  formulated  static  model  for  the 
dynamic  analysis  of  riser  and  developed  an  efficient  algorithm  for  its  numerical  solution. 

O'Brien  and  McNamara  (1988)  developed  a technique  based  on  finite  element 
method  by  the  separation  of  the  rigid  body  motions  from  deformations  of  element  under 
the  condition  of  finite  rotations.  The  method  considered  all  the  nonlinear  effects  of  geo- 
metric changes,  bending-axial,  and  bending-torsional  coupling. 

1.3  Scope  of  Study 

As  mentioned  briefly  earlier,  the  primary  objective  of  this  dissertation  is  to  develop 
a mathematical  model  for  the  dynamic  analysis  of  the  riser  system  with  the  inclusion  of 
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internal  flow  and  nonlinear  effects  due  to  large  structural  displacements  and  solve  it  nu- 
merically for  the  investigation  of  the  effect  of  internal  flow  on  marine  riser  dynamics. 
Iwan-Blevin's  model  is  implemented  as  a vortex-shedding  model  to  examine  the  effect  of 
internal  flow  on  votex-induced  vibration  due  to  inline  current.  A semi-analytical  solution 
using  series  expansion  is  presented  for  obtaining  system  frequencies  and  a conventional 
numerical  method  using  system  matrices  is  employed  to  check  the  results  with  those 
from  semi-analytical  method.  For  the  solution  of  forced  vibration,  Galerkin's  finite  ele- 
ment approximation  is  implemented  and  Iwan-Blevin's  model  is  adapted  to  solve  vortex- 
induced  vibration. 

Chapter  2 formulates  the  mathematical  model  of  the  system.  The  riser  system  is 
modeled  as  a long  tubular  beam  with  the  inclusion  of  internal  flow.  First,  the  simplifica- 
tions and  assumptions  made  in  this  paper  are  described  and  the  equations  of  dynamic 
equilibrium  for  the  differential  element  of  pipe  are  derived  on  the  basis  of  Eulerian  de- 
scription for  the  removal  of  the  complexities  of  the  equations.  The  force  acting  on  the 
interior  wall  of  the  riser,  which  results  from  the  deformation-dependent  accelerations  due 
to  internal  flow,  is  derived  in  three  ways  and  compared  with  each  other.  As  a hydrody- 
namic external  force,  the  normal  force  component  to  the  deformed  riser  is  developed  us- 
ing the  concept  of  relative  motion  between  the  pipe  and  the  surrounding  fluid. 
Combining  all  the  terms,  the  nonlinear  governing  equations  are  derived  with  nonlinear 
boundary  conditions.  The  nonlinear  governing  equations  and  boundary  conditions  are 
approximated  using  the  assumption  of  small  deformation.  Finally,  the  current-vortex 
model  with  the  inclusion  of  internal  flow  is  developed  by  combining  the  approximated 
forms  with  Iwan-Blevin's  model. 

Chapter  3 formulates  the  numerical  model  of  governing  equations  and  boundary 
conditions  using  finite  element  method.  Variational  statement  for  boundary  value  prob- 
lem is  introduced,  that  is,  the  weak  form  of  governing  equation  is  derived.  Sequentially, 
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the  incremental  operator  is  applied  on  the  weak  form  in  order  to  derive  the  incremental 
equilibrium  equation,  and  then  Galerkin's  decomposition  method  is  performed  to  con- 
struct the  approximate  solutions  of  the  problem.  Finally,  incorporation  of  the  usual  cubic 
shape  functions  yields  the  matrix  dynamic  equilibrium  equations  constructed  in  terms  of 
the  unknown  deformations  at  the  nodes.  In  addition,  the  finite  element  formulation  for 
vortex-induced  vibration  is  made  in  the  last  section  using  the  same  procedure  except  ap- 
plying an  incremental  operator. 

In  chapter  4,  the  development  of  a finite  element  code  for  the  solution  of  our  two- 
point  boundary  value  problem  is  discussed.  Some  considerations  of  code  design  are  pres- 
ented and  the  overall  structure  of  the  code  is  outlined. 

The  system  frequencies  and  mode  shapes  of  a marine  riser  with  the  internal  flow 
are  evaluated  in  chapter  5.  This  chapter  consists  of  three  categories;  first,  semi-analytical 
solutions  are  obtained  using  series  expansion  to  the  mathematical  model  and  second,  nu- 
merical solutions  are  calculated  adapting  complex  eigenvalue  method  to  the  matrix  equi- 
librium equation.  Finally,  the  algorithm  for  the  computer  program  is  introduced  and  the 
verification  of  the  first  method  is  inspected  through  the  comparison  of  the  system  fre- 
quencies with  the  second  method.  Further,  the  comparison  of  semi-analytical  solutions 
with  results  from  referenced  papers  is  given. 

Chapter  6 deals  with  the  forced  vibration  of  the  system  derived  in  chapter  3 and  its 
dynamic  behaviors.  The  system  considered  here  consists  of  two;  first,  the  system  simply 
responding  to  the  general  environmental  loads  such  as  wave  or  current  and  second,  the 
system  responding  to  the  vortex  shedding  generated  in  the  transverse  direction  of  inline 
current,  that  is,  vortex-induced  vibration.  Newmark  method  is  used  for  time  step  integra- 
tion and  predictor-corrector  scheme  is  implemented  as  an  iteration  algorithm  for  the  non- 
linear terms.  Using  those  techniques  mentioned  above,  the  vortex-induced  vibration  due 
to  current  and  inline  vibration  due  to  current  or  wave  are  obtained  and  then  their 
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characteristics  are  investigated.  Finally,  the  application  of  top  boundary  condition  to  the 
system  matrices  or  vectors  is  discussed  and  its  effect  on  the  dynamic  responses  of  riser  is 
inspected. 

In  chapter  7,  using  all  the  programs  developed  in  this  study,  the  effects  of  internal 
flow  on  riser  dynamics  are  investigated,  that  is,  the  effects  on  system  frequencies,  vortex- 
induced  vibration,  and  vibrations  due  to  current  or  wave  loading  by  applying  those  pro- 
grams and  numerical  models.  The  investigations  are  carried  out  in  the  way  inspecting 
how  the  internal  flow  affects  the  kinematics  or  stresses  of  the  riser  with  the  various  pa- 
rameters. In  addition,  the  characteristics  of  each  forcing  term  due  to  internal  flow,  name- 
ly, centrifugal  and  coriolis  terms,  and  the  effects  due  to  the  inclusion  of  the  geometric 
and  the  fuidic  nonlinearity  are  discussed  in  this  chapter. 

Finally,  chapter  8 includes  the  conclusions  of  this  study  with  an  investigative  ex- 
planation and  recommendations  for  further  research. 


CHAPTER  2 

MATHEMATICAL  MODEL  OF  THE  SYSTEM 


The  riser  system  can  be  modeled  as  a long  tubular  beam  connecting  the  drilling 
platform  with  the  wellhead  at  the  seabed  and  is  composed  of  rigid  pipes.  Figure  2. 1 de- 
picts the  simplified  configuration  of  the  riser  system.  The  riser  pipes  are  connected  by  the 
riser  connectors.  A ball  joint  is  located  at  the  lower  end  of  the  riser  and  a combination  of 
a ball  and  a slip  joint  at  the  upper  end.  The  ball  joints  alleviate  high  bending  stresses  and 
the  slip  joint  compensates  for  the  heave  motion  of  the  drilling  platform.  A tensioning 
system  is  installed  on  the  drilling  platform  and  applies  a tension  at  the  top  of  the  riser. 
This  tension  provides  part  of  the  support  required  to  keep  the  riser  tight  and  prevent 
buckling  or  collapse.  Additional  supporting  force  is  provided  by  buoyancy  modules  pro- 
perly distributed  along  the  riser.  The  kill  and  choke  lines  are  small  diameter  pipes  with 
high  pressure,  which  run  along  the  riser.  The  lines  are  mounted  directly  on  the  connec- 
tors through  which  they  exert  concentrated  forces  and  moments  of  the  riser. 

The  drill  string  is  inside  the  riser,  which  runs  from  the  drilling  platform  to  the  well. 
The  ring  space  between  the  riser  and  the  drill  string  circulates  the  drilling  mud.  It  exerts 
on  the  riser  static  pressure  force,  coriolis  and  centrifugal  forces  due  to  the  riser's  local 
rotation,  and  vertical  and  torsional  frictional  forces.  And  also,  the  source  of  external 
forces  exerted  on  the  riser  are  the  ocean  currents,  the  surface  and  internal  waves. 

Due  to  the  drifting  of  the  platform  and  the  variation  of  the  current  direction  with 
depth,  the  centerline  of  the  riser  is  in  general  deformed  to  a three-dimensional  curve  with 
the  large  deflections  with  respect  to  its  original  equilibrium  position.  The  mathematical 
model  is  a set  of  consistent  equations  for  the  tubular  beam  with  large  three-dimensional 
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deflections  under  internal  and  external  pressure.  The  extensional  oscillations  of  the  beam 
are  also  included  in  this  formulation. 

2.1  Simplifications  and  Assumptions  of  Model 
The  riser  system  can  be  modeled  as  described  in  Fig.  2.1.  However,  the  inclusion 
of  all  the  terms  mentioned  above  makes  the  analysis  extremely  complicated.  Therefore, 
the  model  derived  in  this  chapter  is  based  on  the  following  simplifications  which  are  usu- 
ally made  in  the  analysis  of  a long  tubular  beam. 


Figure  2.1:  Simplified  Configuration  of  the  Riser  System 
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1 . The  riser  can  be  modeled  as  a beam  rather  than  a shell  because  the  ratio  of  its 
diameter  to  the  riser  length  is  small.  If  detailed  analysis  is  necessary  locally,  the  riser  can 
be  locally  modeled  as  a shell.  However,  not  only  our  knowledge  of  the  distribution  of 
hydrodynamic  loads  around  cylinder  is  very  limited  but  also  detailed  analysis  is  also 
complicated. 

2.  The  drill  string  occasionally  comes  in  contact  with  the  riser  through  the  tool 
joints  and  exerts  concentrated  forces  and  moments  on  the  riser.  The  presence  of  the  drill 
string  is  ignored.  This  will  underestimate  the  weight  of  the  riser. 

3.  The  frictional  force  exerted  on  the  riser  due  to  the  mud  flow  is  not  taken  into 
account.  Consequently,  the  plug-flow  model  with  no  radial  variation  of  velocity  is 
utilized. 

4.  The  variation  of  temperature  of  water  and  drilling  mud  with  depth,  which  may 
induce  stresses,  is  neglected. 

In  addition  to  the  foregoing  simplifications  of  the  riser  system,  the  following  as- 
sumptions are  made. 

1.  The  riser  material  is  isotropic,  homogeneous,  elastic,  and  linear,  i.e.,  it  obeys 
Hook's  law. 

2.  Plane  sections  remain  plane  after  deformation,  i.e.,  warping  is  neglected. 

3.  Cross  sections  remain  normal  to  the  neutral  axis  after  deformation,  i.e.,  shear 
deformation  is  neglected.  Thus,  the  riser  system  can  be  modeled  as  a Bemoulli-Euler 
type  of  beam  and  not  as  a Timoshenko  type. 

4.  The  initial  undeformed  configuration  of  the  riser  is  straight  even  though  it  is  ver- 
tically inclined,  i.e.  initial  imperfections  are  not  present.  Therefore,  both  initial  curva- 
tures and  geometric  torsion  are  zero,  and  adjacent  cross  sections  of  the  riser  are  parallel 
to  those  in  the  undeformed  state. 
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Based  on  these  simplifications  and  assumptions,  the  model  for  a three-dimensional 
dynamic  behavior  of  the  riser  is  established.  It  is  a set  of  consistent  equations  for  beam 
with  large  deflections  and  consider  the  extensional  oscillations  of  the  riser.  The  hydrody- 
namic load  due  to  wave  and  current  is  integrated  over  the  external  surface  of  the  riser  and 
the  fluid  force  due  to  internal  flow  is  exerted  over  the  internal  surface  of  the  riser.  The 
boundary  conditions  can  be  changed  by  slightly  modifying  the  model  and  different  initial 
conditions  can  be  imposed. 

2.2  Equations  of  the  Pipe 

This  section  derives  the  equations  of  dynamic  equilibrium  for  a differential  element 
of  the  pipe.  The  centerline  of  the  riser  is,  in  general,  a three  dimensional  curve  and  Euler- 
ian  description  and  vector  representations  are  applied  for  the  removal  of  the  complexities 
of  the  equations. 

The  main  aspects  of  the  curved  geometry  of  the  system  are  depicted  in  Figure  2.2. 
Three  coordinates  systems  are  defined,  namely 

AAA 

• ( ‘,j>  k),  an  orthonormal  global  inertial  system  of  coordinates  with  origin  at  the 
lower  ball  joint  of  the  riser. 

A A A 

• ( n,  b,  t ),  an  orthonormal  triad  passing  through  the  centroid  of  the  cross  section 
of  the  deformed  configuration. 

AAA 

• (N,  B,  T),  an  nonorthonormal  triad  passing  through  the  centroid  of  the  cross  sec- 
tion of  the  deformed  configuration. 

where  "A"  indicates  a unit  vector  and  hereafter  boldface  represents  general  vectors. 

Since  we  assume  no  shear  deformation,  the  third  nonorthonormal  triad  becomes 
consistent  with  the  second  orthonormal  triad.  The  initial  undeformed  configuration  of  the 

A 

riser  is  vertically  straight,  i.e.,  its  direction  accords  with  k.  The  deformed  space  curve  is 

AAA 

specified  by  giving  its  position  vector  r from  the  origin  of  space  fixed  system  ( i,j,  &)  to 
a point  on  the  deformed  riser  centerline  ( see  Fig.  2.2  ),  as  a function  of  deformed 
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arclength  s . At  any  point  on  the  deformed  curve  the  unit  tangent  vector  t,  the  unit  nor- 
mal vector  h,  and  the  unit  binomial  vector  b are  defined  by 

A . AA/  A/VA 

t = r , n - r /k  , b = t x n ( 2.1  ) 


where  prime  denotes  differentiation  with  respect  to  s , and  the  curvature  k,  from  the  well 
known  results  about  the  differential  geometry  of  curves  in  space  (Love,  A.,  1927),  is  giv- 
en by 


K2  = 


„//  „// 


or 


- r 


„/// 


(2.2) 


where  " . " denotes  inner  product  of  vectors  ( dot  product ). 


Figure  2.2:  Position  Vector  and  Principal  Direction  at  a Point  on  the  Riser  Centerline 
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The  total  bending  moment,  as  shown  in  Fig.  2.2,  acts  in  the  binomial  direction  and 
is  proportional  to  the  bending  rigidity  of  the  cross  section  El,  and  the  sign  of  curvature  is 
fixed  by  convention  since  change  of  sign  merely  reverses  the  direction  of  h and  b as  seen 
from  Eq.  (2.1).  For  three  dimensional  curves,  it  is  convenient  to  choose  the  positive  sign 
so  that  n is  in  a fixed  direction  normal  to  the  plane  and  (h  , b , t)  form  a right-  handed 
orthonormal  triad.  Hereafter,  the  space  curve  will  be  identified  with  the  central  axis  of 
the  riser  in  the  deformed  state,  with  the  position  vector,  as  shown  in  Fig.  2.2,  represented 
as  follows. 

AAA  A 

r = r„  + sk  = Xi  i + x2j  + ( s + X3  ) k ( 2.3  ) 

where  r0  is  the  deformation  vector  at  any  point  on  the  riser  in  the  undeformed  state. 

In  the  classical  theory  of  rods,  the  internal  state  of  stress  at  a point  on  the  rod  is  ful- 
ly characterized  by  the  resultant  force  F and  the  couple  moment  M acting  on  the  central 
axis.  Conservation  of  linear  momentum  leads  to  the  equation  of  motion 

F'  + Fe  + Fi  = mrr  ( 2.4  ) 

where  Fj  and  FE  are  the  internal  fluid  force  and  the  applied  external  hydrodynamic 
force  per  unit  length,  respectively  (see  Fig  2.2)  and  mr  is  the  mass  of  the  riser  per  unit 
length,  and  dot  denotes  differentiation  with  respect  to  time.  Since  the  centerline  of  the 
riser  is  deformed  to  a three  dimensional  curve,  the  conservation  of  angular  momentum 
must  include  the  rotatory  inertia  even  if  the  shear  deformation  is  negligible.  Thus  the  mo- 
ment conservation  equation  becomes 

M/  + /xF  + m = [Jr]  {to}  ( 2.5  ) 

where  m is  the  distributed  couple  vector  induced  by  the  asymmetric  flow  due  to  vortex 

shedding.  It  is  usually  negligible  if  we  ignore  the  coupling  of  bending  and  torsion  is 
weak.  The  mass  moment  of  inertia  [Jr]  is  given  by 
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[Jr]  = 


Jn  0 0 

0 J22  0 

0 0 J33 


and  to  is  the  angular  acceleration  and  is  obtained  from  the  time  differentiation  of  the  an- 
gular velocity.  The  angular  velocity  to  is  computed  using  the  appropriate  transformation 

AAA  A A 

matrices  that  carry  the  initial  system  ( / , j , k)  to  the  final  triad  ( h , b , 1 ).  Transforma- 
tion from  the  initial  system  to  the  final  triad  takes  place  through  three  translations 
( , *2 , *3  ) and  three  rotations  ( 9i , 62 , 63  ) in  the  following  sequence: 

A AAA 

• rotation  03  about  the  k axis  to  get  the  intermediate  system  ( / 1 , j 1 , k\  = k) 

A A 

• rotation  02  about  the  new  axis  j \ to  get  the  second  intermediate  system  ( /' 2 , 

A A A 

j 2 =j  1 , k2  ) 

A A A 

• rotation  0i  about  the  new  axis  17  to  get  the  final  triad  ( h , b , t ) 

As  a results,  the  angular  velocity  m is  given  by 

tu  = 63^1  + 02^2  + 0i  h ( 2.6  ) 


To  perform  the  above  rotations,  we  have  to  apply  sequentially  the  following  rotation  ma- 


trices 


COS  03 

sin03 

0 

[R>]  = 

-sin  03 

COS  03 

0 

0 

0 

1 

COS  02 

0 

-sin  02 

[R2]  = 

0 

1 

0 

sin  02 

0 

COS  02 

1 

0 

0 

[R3]  = 

0 

COS  01 

sin  0i 

0 -sin  0i  cos©i 


(2.7) 


(2.8) 


(2.9) 
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Using  Equations  (2.7),  (2.8)  and  (2.9),  Eq.  (2.6)  can  be  presented  in  terms  of  the  fixed 
inertial  system 


tn  = ( 0icos02cos03  - 02sin03  J / + ( 0icos02sin03  + 02cos03  1 j 


+ 103-  0isin02  j k 


(2.10) 


The  angular  acceleration  to  can  then  be  obtained  by  time  differentiation  of  Eq.  (2.10). 

According  to  the  classical  Bemoulli-Euler  theory  of  elastic  rods  with  equal  princi- 
pal stiffness  (Love,  A.E.H.,  1927),  the  bending  component  of  the  stress  couple  is  propor- 

A 

tional  to  the  curvature  k and  is  directed  along  the  binomial  vector  b.  In  addition,  the 

torsional  component  H of  the  stress  couple  is  directed  to  the  tangent  vector  It . Thus,  the 
constitutive  equation  for  the  stress  couple  can  be  written  as 

M = EIk£  + H / (2.11  ) 


Upon  substitution  of  Eq.  (2.11)  into  Eq.  (2.5)  after  elementary  vector  manipula- 
tions using  Eq.  (2.1),  we  have 


- HkA  + F 


+ H'i  + m = [Jr]  {to} 


(2.12) 


The  dot  and  cross  product  of  t with  Eq  (2.12)  using  Eq.  (2.2)  yield  respectively 


H'  + t-  m = It  ■ ([JR]  {to})  ( 2.13  ) 

F = - (El r//)/  + (Te  - EIk2) r'  + r'  x m + Hr'  x r"  - r'  x [Jr]  {to}  ( 2.14  ) 


where  Te , the  effective  riser  tension,  is  defined  by 


Te  = t-¥ 


(2.15) 
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Alternatively  the  effective  riser  tension,  which  is  the  tangential  component  of  the  inter- 
nal force,  may  be  given  approximately  by  (Dareing,  1976) 

Te  * T + pwgnD20  (Sw  - z)/4  - pmgnD?  (Sm  - z)/ 4 ( 2.16  ) 

where  T = actual  tension  in  the  riser,  Sw  = ordinate  of  the  free  surface  of  the  water,  Sm  = 
ordinate  of  the  free  surface  of  mud,  pw  and  p,„  are  the  densities  of  water  and  mud  re- 
spectively, D0  and  Di  = external  and  internal  diameter  of  the  riser  and  z is  measure  from 
the  lower  ball  joint,  i.e.,  z = s + x3. 

The  actual  tension  T also  satisfies  the  constitutive  relation: 


T = EAe,  (2.17) 

where  EA  is  the  stretching  rigidity  and  e,  is  the  strain  of  the  riser  centerline  in  the  tan- 
gential direction  defined  as 


T dr_  dx_  ~\  2 

Jt  L ds  ds  J 


where  s is  related  to  s,  by 
ds  i 


— = 1 + S,  or  — = 


ds 
ds\ 


1 + 8f 


(2.18) 


(2.19) 


For  an  inextensible  riser,  s,  = 0 and  the  deformed  arclength  s,  becomes  the  undeformed 
arclength  s from  Eq.  (2.19).  However,  Eq.  (2.16)  and  Eq.  (2.17)  have  to  be  kept  to  com- 
pute £/  for  the  extensible  riser. 

Equation  (2.13)  indicates  the  torsional  relation.  If  no  distributed  torsional  couple  is 
applied  to  the  riser  and  a rotatory  inertia  is  neglected,  as  shown  in  Eq.  (2.13),  H is  inde- 
pendent of  s,  just  as  in  the  static  theory  of  rods  with  equal  principal  stiffnesses.  If  a dis- 
tributed torsional  couple  acts  or  a rotatory  couple  is  included,  then  H can  be  determined 
by  Eq.  (2.13)  and  the  prescribed  value  of  H at  one  end. 
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Substituting  Eq.  (2.14)  into  Eq.  (2.4),  we  can  obtain  the  vector  equation  of  motion 
for  the  riser 

- (El r//)//  + [(Te  - EIk2)  r7]7  + (r7  x m)7  + (Hr7  x r77)7  - (r7  x [Jr]  {tn})7 

+ Fe  + Fi  = mr  r ( 2.20  ) 

In  this  equation,  the  structural  damping  term  is  not  considered.  The  effective  tension  Te 
is  unknown  and  r must  also  satisfy  the  extensibility  or  the  inextensibility  condition  as 
discussed  earlier. 

2.3  Force  due  to  Internal  Flow 

In  this  section,  the  force  acting  on  the  internal  wall  of  the  riser  is  derived  in  three 
ways:  by  obtaining  the  acceleration  of  a fluid  particle,  by  using  the  concept  of  Hamilton’s 
principle,  and  by  evaluating  the  forcing  components.  In  all  three  ways,  no  small-scale 
motions  such  as  turbulence  or  secondary  flow  are  assumed  to  be  absent.  And  also,  the 
plug-flow  model  with  no  radial  variation  of  velocity  is  utilized  as  a fluid  model  for  the 
internal  flow  as  discussed  in  section  2. 1 . 

First,  we  obtain  the  acceleration  of  a fluid  particle  by  differentiating  the  fluid  ve- 
locity Vj, 

cN  i dVi  ds  dV\ 

Al  = ~dT  = ~~  + ~ (2'21) 

where 

dx\  a dx2  a dx-x  * A 

Vl  = + ~j  + ~k  + VH  = r + Vrr7  ( 2.22 ) 

Performing  the  operations  in  Eq.  (2.21)  with  s\  &s  and  ds/dt  = Vi,  we  obtain  the  fluid 
acceleration  in  the  deformed  direction 

Ai  = (r  + VirO'Vi  + Vi 

= r + Vir7  + 2 Vi  r7  + ViV(r7  + V\r" 


( 2.23  ) 
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Assuming  steady  flow  and  no  convective  acceleration  along  the  riser  curve,  the  fluid 
force  acting  on  the  internal  wall  of  riser  becomes 

Fi  = - tfi/Ai  = - mf(  r + 2Vir/  + Vfr"  ) (2.24) 

Second,  we  apply  the  concept  of  Hamilton's  principle  just  considering  the  finite 
part  comprising  the  pipe  and  the  enclosed  volume  of  fluid.  Hamilton's  principle  in  our 
problem  states  that 

5 


f 


(Tr  + Tf  - Vr  — V/) dt  — 0 


( 2.25  ) 


where,  Tr  and  Vr  are  the  kinetic  and  potential  energies  associated  with  the  tube,  and  Tf 
and  Vf  are  the  corresponding  quantities  for  the  fluid.  Using  the  velocity  of  the  internal 
fluid  obtained  from  neglecting  the  stretching  strain  in  Eq.  (2.22),  we  have  kinetic  energy 
of  the  enclosed  volume  of  fluid 


dx\  dxx 

Vi  + ( dt  +Vl  + ( 


dx2  dx2 
~+  Vi  — )2 


dt 


ds\ 


dx3 


dx3 


+ ( a/'+Vl  ds]  + l )2  l*1 


( 2.26 ) 


and  the  potential  energy  of  the  fluid  is  zero  because  the  fluid  is  assumed  to  be  incom- 
pressible, i.e., 

Vf=  0 ( 2.27  ) 

Performing  the  variation  after  the  substitution  of  Eq.  (2.26)  and  (2.27)  into  (2.25),  we 
obtain  the  following  integration 

j^[>/{(JCi  + VIx/1)(8xi  + VjSx' ) + (x2  + Vi*2  )(8x2  + ViSxj) 

+(jc3  + V1X3  + 1 )(8x3  + ViSxj  ) } + 87V  - bVr\dsi  dt  ( 2.28  ) 
5jc;  = (5x/) , 8x-  = —.(Sxj) 


Since 


(/-  =1,2,3) 
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each  terms  may  be  integrated  by  parts  so  as  to  eliminate  the  various  derivatives  of  5x,-. 
When  this  is  done,  there  is  obtained 

T f [ E ( nt/Xi  + 2/w/Vix'  + + 5 Tr  - bVr]ds\  dt  - 0 ( 2.29  ) 

J( i •>/  f=i 

where  all  the  integrated  terms  have  disappeared  because  of  the  boundary  conditions. 
From  the  concept  of  the  resulting  Euler-Lagrange  equations,  we  can  recognize  that  the 
expression  in  round  bracket  represents  the  forcing  components  of  a fluid  particle  inside 
the  pipe  due  to  internal  flow.  Thus,  the  fluid  force  acting  on  the  internal  wall  of  the  pipe 
can  be  written  in  vector  form 

Fi  = - /»/( ir  + 2 Vi  r/  + Vj  r"  ) (2.30) 

We  can  easily  see  that  Equations  (2.24)  and  (2.30)  are  identical.  The  first  term  on  the 
right  represents  the  inertia  force  associated  with  the  riser  acceleration.  The  second  term  is 
the  inertia  force  associated  with  the  coriolis  acceleration  which  arise  because  the  fluid  is 
flowing  with  velocity  V,  relative  to  the  riser,  while  the  riser  itself  has  an  angular  velocity 
at  any  point  along  its  length.  The  last  term  represents  the  inertia  force  associated  with  the 
change  in  direction  of  the  flow  velocity,  owing  to  the  curvature  of  the  riser. 

Finally,  we  can  evaluate  all  the  forcing  components  in  a physical  sense.  There  are 
three  kinds  of  inertia  forces  by  the  drilling  fluid:  the  force  due  to  linear  acceleration  of 
the  riser,  the  coriolis  force,  and  the  centrifugal  force.  From  the  general  knowledge  of 
those  components  (Hibbler,  R.C.,  1986),  they  are  given  by  Eq.  (2.31)  to  (2.33), 
respectively. 


Fi  ia  = mfr 

(2.31) 

Fico  = 2/w/nj  x Vi 

( 2.32  ) 

F ice  = mc  x (m  x h! k) 

( 2.33  ) 

A A 

where  tnc  = Vi  k b and  Vi  = Vi  t , and  also  tdc  can  be  calculated  from  Eq.  (2.10). 
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2.4  Hydrodynamic  Forces 

The  relative  motion  between  the  pipe  and  the  surrounding  fluid  produces  hydrody- 
namic force  composed  of  inertia,  drag,  and  frictional  forces.  The  resultant  total  force  dis- 
tribution Fe  along  riser  length  is  decomposed  into  a normal  force  component,  Fn,  and  a 

tangential  component,  Ft..  To  perform  this  decomposition,  the  relative  velocity  and  ac- 
celeration is  resolved  into  components  perpendicular  and  parallel  to  the  deformed  riser 

axis  as  shown  in  Figure  2.3.  The  relative  velocity  vector  between  the  pipe  and  the  sur- 
rounding fluid  is  given  by 

Vr  = Vp-Vc-Vw  - r -Vc-Vw  (2.34) 


where  Vp  = pipe  velocity  vector,  Vc=  steady  current  velocity  vector,  and  Vw  = time  de- 
pendent wave  velocity  vector.  It  should  be  noted  that,  current  and  wave  velocity  are  func- 
tions of  vertical  position.  Figure  2.3  depicts  the  velocity  vectors  acting  on  a segment  of 
pipe  and  explains  the  definition  of  the  relative  velocity  components.  Noting  that  the  cur- 
rent velocity  is  independent  of  time,  the  relative  acceleration  of  the  pipe  with  respect  to 
the  fluid  is  given  by 


Vr  = Vr  - Vw  = r - Vw 


(2.35) 


Pipe  velocity: 

Fluid  velocity  (current  + waves): 
Relative  velocity: 

Pipe  tangent  vector: 

Tangential  component  of  Vr: 
Normal  component  of  Vr: 


Vp  = r 
VF 


Vr  = r - VF 


t 


VrT  = ( V R ■ t)t 
Vrn  = Vr  - Vrt 


Figure  2.3:  Velocity  Vectors  Acting  on  a Segment  of  Pipe  and  Definition  of  the  Rela- 
tive Velocity  Components 
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and  its  tangential  and  normal  components  are  given  by 

Vrt  = (VR  /)/  and  Vrn  = Vr  — Vrt  ( 2.36 ) 

Also,  we  need  the  tangential  and  normal  components  of  the  wave-induced  water  particle 
acceleration  to  calculate  wave  induced  force  and  they  are  given  by 

V wx  = (Vw-o*  and  V wn  = Vw  - V wt  (2.37) 

The  normal  component  of  the  relative  velocity,  the  current  velocity,  and  the  wave  veloc- 
ity vectors  are  expressed  as  a component  form  in  xi,  X2,  x3  directions. 

AAA 

Vrn  = Vrni  i + Vrn27  + V RN3  k 

A A A 

Vc  = uc/  + vcy  + wcA: 

A A A 

Vw  = uw  / + vwy  + ww  k ( 2.38  ) 

where, 

Vrni  = Xl-Uc-Uw  - [x\  (Xi-Uc-Uw)  + Xj^-Vc-Vw) 

+ (x3  + 1/(1  +e,))(x3-wc-ww)]x/1 

V RN2  = *2-Vc-Vw  “ (il-Uc-Uw)  + Xj  (X2-Vc-Vw  ) 

+ (Xj  + 1/(1  +S,))(x3-Wc-Ww)]X2 

Vrn3  = X3-Wc-Ww  - [Xj  (Xi-Uc-Uw)  +X2(X2-Vc-Vw) 

+ (X3  + 1/(1  +S,))(X3-W0-WW)](X3  + 1/(1  + 8,)) 

(2.39) 

and  the  normal  component  of  the  relative  acceleration  and  the  wave  velocity  vectors  as  a 
in  directional  components 

A . A A 

Vrn  = Vrni  i + Vrn27  + Vrn3  k 


Vw  = uw/  + vwy  + wwA: 


( 2.40 ) 
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where, 

[x\  (*i-uw)  +x'2(x2-Vw)+(x3  + 1/(1  + e,))(x3-ww)]x/1 
[x'  (Xi-Uw)  + X2  (*2-Vw)+(*3  + 1/(1  + e,))(x3-Ww)]x2 
[x[  (Xi-Uw)  + X2(x2-Vw) 

+ (X3  + 1/(1  +8,))(X3-WW)](X3  + 1/(1  +8,))  (2.41  ) 

Finally,  the  normal  component  of  the  water  particle  acceleration  of  wave  is  given  by 

VwN  = V WNl  i + VwN2  j + VwN3  k ( 2.42  ) 

where, 

V WNl  = Uw  - [XjUw  +X2VW  + (x3  + l/(l  +S,))WW  ] Xj 
V WN2  = Vw  - [XjUw  +X2VW  +(x3  + 1/(1  +8/))  Ww  ] X2 

VwN3  = Ww  - [x',  Uw  + X2  Vw  +(X3  + 1/(1  +8,))WW](X3  + 1/(1+  8,))  ( 2.43  ) 

the  wave  or  particle  velocity  and  acceleration  are  obtained  from  the  velocity  potential 
approach. 

In  the  present  case,  the  tangential  force  component  by  skin  friction  is  neglected, 
and  only  the  normal  force  component  is  retained  by 

Fe  ~ Fn  = - pw  (tlDq/4)Ca  Vrn 

+ pw  (7iDq/4)  Vwn  - PwDo  Cd  |Vrn|Vrn/2  ( 2.44  ) 

where  |Vrn|  = [V^  + V^,2  + VRN3]1/2  and  CA  = added  mass  coefficient. 

The  first  term  of  Eq.  (2.44)  is  an  inertia  term  representing  the  force  that  is  required 
to  accelerate  the  pipe  with  respect  to  the  surrounding  water.  The  second  term  is  the  wave 
induced  force.  This  force  is  produced  by  the  local  pressure  gradient  that  accompanies  the 
normal  component  of  water  particle  acceleration.  The  last  term  is  the  drag  force  that  is 


Vrni  = Xi-Uw  - 
VrN2  = X2-vw  - 

Vrn3  = X3— Ww  - 
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proportional  to  the  square  of  the  normal  velocity  component  and  is  formed  by  the  separa- 
tion of  flow. 

2.5  Equations  of  Motion  and  Boundary  Conditions 
2.5.1  Nonlinear  Governing  Equations 

Eq.  (2.20)  is  the  bending  equations  of  motion  of  the  pipe.  In  this  equation,  H is  the 
time  dependent  torque  due  to  variation  of  the  relative  orientation  of  the  supporting  plat- 
form and  imperfections  of  the  tensioning  system,  and  m is  the  distributed  couple  vector 
induced  by  the  asymmetry  of  the  flow  due  to  vortex  shedding.  But  m is  usually  negli- 
gible, leading  to  the  uncoupling  of  the  bending  and  torsion.  Furthermore,  rotatory  inertia 
term  can  be  neglected  because  the  riser  response  is  usually  dominated  by  low  frequency 
motion.  Thus,  one  deduces  from  Eq.  (2.13)  that  H is  independent  of  the  pipe  length,  that 
is,  it  has  constant  value  along  the  riser  arclength.  Eq.  (2.20)  then  becomes 

-(El r")"  + [(Te  - EIk2)  r/]/  + H(r'  x r")'  + FE  + Fi  = mr  r ( 2.45  ) 

Substituting  Eq.  (2.30)  and  (2.44)  into  the  above  equation 

mr r + (Elr")"  - [(Te  - ElK2)r']'  - H(r'  x r")' 

= - nif(  r + 2 Vif'  + Vj  r"  ) 

- pw  (7tZ)o/4)  Ca  Vrn  + p w (7tDo/4)  Vwn  - PwDo  Cd  |Vrn|  Vrn/2 

(2.46) 

Manipulating  and  rewriting  the  last  three  terms  of  right  hand  side  after  the  substitution  of 
Eq.  (2.36)  and  (2.37)  into  the  last  three  terms  of  (2.46) 

- ^PwCdD0|Vrn|Vrn  - Pm-AoCai1  + P»-A0(1+Ca)Vw 

+ Pw>A0Ca(^  ■ — pw Ao(l+C a)(^  ■ V w)^  ( 2.47  ) 


where  A0  = tiD20/4 
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Substituting  Eq.  (2.47)  into  (2.46),  and  letting 

q = - ^PwCdD0|Vrn|Vrn  + PwA0(1+Ca)Vw 
+ PwA0Ca(?  • r)^  - pwA0(1+Ca)(^  • Vw)/ 

and  mt  - mr  + nif  + pwA0CA  ( 2.48  ) 

we  have  nonlinear  governing  equations  (N. G.E.'s). 

mtr  + ImfVir'  + m/V^r"  + (Elr")^ 

- [(Te  - ElK2)r/]/  - H(r/  x r")'  = q ( 2.49  ) 

The  effective  tension  needed  in  the  above  equation  can  be  computed  in  two  ways: 
first  by  using  Eq.  (2.15)  and  unknown  internal  force  F found  from  the  local  equilibrium 
equation  which  will  be  derived  in  next  chapter  and  second  by  using  Eqs.  (2.16)  and 
(2.17).  Actually,  the  first  method  should  be  used  to  evaluate  Te  and  the  second  to  com- 
pute £<  for  extensible  risers.  For  inextensible  risers,  e,  = 0 must  be  used  instead  of  Eq. 

(2.17). 

Utilization  of  the  kinematic  constraint  on  the  unit  tangent  vector,  |/|  = |r'|  = 1, 
leads  to  the  reduction  of  number  of  the  unknowns.  The  following  kinematic  relations  are 
produced  from  this  constraint. 

*3  - (l-x\2 -x'22y/2  - 1/(1 +E/)  (2.50) 

x3(vi)  = x3(0)  + Jo1  Xjdsi  (2.51) 

Consequently,  the  vertical  deflections  can  be  computed  in  terms  of  lateral  one.  The  kine- 
matic relations  are  used  as  predictors  to  remove  the  vertical  unknowns  from  the  govern- 
ing equation.  They  are  used  again  as  correctors,  after  the  solution  of  the  reduced 
equations  has  been  obtained.  This  will  reduce  the  iteration  or  time  steps. 
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Also,  the  derivative  of  the  position  vector  with  respect  to  time  and  deformed  ar- 
clength  can  be  represented  by 


2.5.2  Boundary  Conditions 

For  the  completion  of  the  mathematical  model,  the  boundary  conditions  have  to  be 
specified.  The  boundary  conditions  for  the  dynamic  model  of  riser  are  usually  time  de- 
pendent because  of  the  motions  of  the  supporting  platform.  The  riser  support  can  be 
modeled  by  substituting  the  linear  translational  or  rotational  springs  providing  the  restor- 
ing boundary  force  and  moments.  Typically  the  displacement  or  the  force  vector  and  the 
unit  tangent  vector,  or  the  bending  or  torsional  moment  has  to  be  specified  at  each  riser 
end.  Figure  2.4  describes  the  free  body  diagram  about  the  equilibrium  of  forces  for  the 
differential  element  at  the  top  of  the  riser.  From  the  free  body  diagram,  equilibrium  of 
forces  at  the  top  yields: 


It  = r'  = x[  i + x'2j  + (*3  + 1/(1  + et))k 


AAA 


r = Jci /'  + *2  j + X3  k 


(2.52) 


Fiu  = (TF)  t ■ i - K\  xi„ 


( 2.53  ) 


F2u  = (TF)  t j - K2x2u 


(2.54) 


F3u  = (TF)  t k-  K2x3u  + TTR 


(2.55) 


where 


TF  » pwgnD20(Sw  - z)/4  - pmgnD:j(Sm  - z)/ 4 


(2.56) 
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the  subscript  u indicates  the  upper  end  of  the  riser,  TTR  is  the  tension  applied  at  top  of 
the  riser  by  the  tensioning  system  and  K,,  K2,  K3  are  spring  constants  supplied  by  the  re- 
storing boundary  force.  Equations  (2.53)  - (2.55)  show  that  the  fluidic  tension  at  the  top 
of  the  riser  always  acts  in  the  tangential  direction.  Both  TF  and  its  direction  are  deforma- 
tion dependent.  Further,  the  internal  force  F is  also  deformation  dependent.  The  nonlin- 
earity of  these  boundary  conditions  disappear  in  linear  model  because  the  three  direction 
cosines  in  (2.53)  - (2.55)  become  0,0,  and  1 respectively. 

The  fixed-hinge  support  are  applied  at  the  lower  end  of  the  riser  as  boundary 
conditions.  The  reaction  force  vector  or  zero  displacement  vector  and  zero  moment  or  the 
unit  tangent  vector  are  given,  and  the  torsional  moment  at  the  lower  or  upper  end. 


Figure  2.4:  Free  Body  Diagram  of  the  Equilibrium  of  Forces  and  Moments  for  the  Dif- 
ferential Element  at  the  Top  of  the  Riser 
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2.5.3  Approximate  Governing  Equations  and  B.C.s 

In  two  previous  sections,  the  geometric  nonlinearities  resulted  from  large  deflec- 
tion and  rotation  have  been  considered  in  the  derivation  of  the  governing  equations  and 
boundary  conditions.  Eulerian  description  has  been  implemented  to  remove  the  complex- 
ity of  the  equations  and  the  differentiation  with  respect  to  the  deformed  arclength  have 
been  operated  in  their  expressions.  However,  it  is  no  longer  necessary  for  the  analysis  of 
the  riser  in  an  intermediate  water  depth  and  we  can  utilize  an  assumption  of  the  small 
deflection  beam  theory. 

In  the  riser  problem  in  intermediate  water  depth,  the  deformed  arclength  Sj  can  be 
approximately  identified  with  the  undeformed  arclength  s for  inextensible  riser  and  the 
nonlinear  vector  equation  (2.20)  can  be  used  except  the  fact  that  prime  denotes  differenti- 
ation with  respect  to  s (Nordgren,  R.P.,  1974).  Further,  if  we  notify  the  approximation 
used  to  obtain  the  fluid  acceleration  due  to  internal  flow  (between  Eqs.  (2.22)  and 
(2.23)),  we  can  have  the  following  equation: 


where  prime  denotes  the  differentiation  with  respect  to  the  undeformed  arclength  s. 
For  the  extensible  riser,  the  actual  tension  is  given  by  Hook's  law  as 


For  nearly  vertical  riser,  we  introduce  a rectangular  Cartesian  coordinate  system 
with  the  x3  axis  vertically  straight  upward.  And  then,  with  small  deflection  and  rotation, 

i.e.,  x\ , x'2  ((  1 and  s + x3  » s,  the  curvature  and  the  strain  along  centerline  are  given  by 


mtr  + 2w/Vif/  + zw/Vf  r"  + (Elr")^ 

- [(Te  - EIk2)t/]/  - H(r'  x r")'  = q 


(2.57) 


T = EAs,,  8/  = (r'V)1/2  - 1 


(2.58) 


■3/2  \2  1 1/2 


(2.59) 
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&t  — x$  + (x^2  + x*2  )/2  (2.60) 

In  the  above  equation,  the  nonlinear  terms  in  bracket  have  been  retained  in  order  to  in- 
clude the  effect  of  lateral  deflection  on  tensions. 

Restricting  attention  to  risers  whose  deformation  lies  wholly  in  a single  plane  and 
dropping  nonlinear  terms,  the  vector  equation  ( 2.57  ) reduces  to 

m,x  1 + 2m /V \x\  + m/V \x'[  + (EIxf)//  - (Tex')'  + Hx"'  = qi  (2.61  ) 

m,x 2 + 2/W/V1X2  + m/V  1X2  + (EIx"  )"  - (TeXj)'  - Hxf'  = q2  (2.62) 

m,x3  - Ti  = q3  ( 2.63  ) 

Most  riser  problems  in  intermediate  water  depth,  longitudinal  vibration  is  unimpor- 
tant and  we  may  neglect  the  longitudinal  inertia  term  in  Eq.  (2.63).  Further,  it  is  conve- 
nient to  include  the  hydrostatic  effects  of  internal  and  external  fluid  pressures  by  defining 
effective  weight  per  unit  length  w and  effective  tension  Te  as 

w — wr  + Y/A/  - y0A0  ( 2.64  ) 

Te  = T ~PiAj  + P0A0  ( 2.65  ) 

where  wr  =riser  weight  per  unit  length,  y,- , y0  = specific  weight  of  internal  and  external 
fluid,  pi  ,p0  = internal  and  external  static  pressure  of  riser,  and  A,- , A„  = internal  and 
external  area  of  riser. 

Eq.  (2.64)  leads  to  q3  = - w.  Eq.  (2.63)  can  then  be  integrated  to  give 


Te  = TTR  - j0  wds 

(2.66) 

Te  = TTB  + ^ wi/ds 

(2.67) 

where  TTR  and  TTB  are,  respectively,  the  top  tension  and  the  bottom  reaction  tension. 
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Thus,  the  effective  tension  is  independent  of  the  horizontal  deflections  and  then, 
with  no  torsion,  Eqs.  (2.61)  and  (2.62)  becomes 

m,x i + 2mfV\x[  + m/V \x'{  + (El x'{)"  - wx'  - Tex'{  = qi  (2.68) 

mtx2  + 2m/Vix'2  + mfV \x2  + (EIx")"  - wxj  - Tex"  = q2  (2.69) 

Equations  (2.68)  and  (2.69)  represent  the  linear  model  and  can  be  also  derived  by 
applying  Hamilton's  principle  (Chen,  B.C.M.,  1992).  Detailed  are  given  in  Appendix  A. 

The  boundary  conditions  associated  with  linear  model  can  be  found  by  modifying 
the  direction  cosines  of  Eqs  (2.53),  (2.54)  and  (2.55)  into  0,0,  and  1 and  are  given  as 
following: 

Fiu  = — K\  X\u 

F 2u  = - K2x2u 

F3u  = TF  - K3x 3„  + TTR  ( 2.70  ) 

2.5.4  Current- Vortex  Model 

In  this  section,  the  mathematical  model  of  vortex-induced  vibration  with  internal 
flow  is  presented  using  the  two  dimensional,  coupled  wake  oscillator  model  proposed  by 
Iwan  and  Blevins  (1974)  for  the  Reynolds  number  of  103to  105.  There  are  many  vortex 
excitation  model  such  as  harmonic  model,  wake  oscillator  model,  etc.  Of  those  models, 
Iwan-Blevin's  model  has  a advantage  as  a numerical  model.  This  model  has  a nature  of 
self-excited  vortex  shedding  which  means  that  the  fluid  behavior  may  be  modeled  by  a 
simple,  nonlinear,  self-excited  oscillator.  Attention  is  confined  to  plane  vibrations 
(x2  = 0)  transverse  to  a steady  current  in  direction  x2 . In-line  vibration  is  treated  in  the 
usual  way  and  no  coupling  is  considered  in  this  model  (see  Figure  2.5). 

First,  a so-called  hidden  flow  variable  w is  defined  as  a fluid  motion  transverse  to  a 
steady  current  and  it,  in  turn,  produces  transverse  force  to  cylinder.  This  transverse  force, 
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which  induces  the  motion  of  the  cylinder  is  a forcing  component  in  the  oscillator  equa- 
tion. The  net  force  per  unit  length  on  cylinder  is  evaluated  using  the  momentum  equation 
for  the  control  volume  containing  cylinder  as  shown  in  Figure  2.5  and  is  given  by 

qi  = a4pwD0uc(w  - xx)  ( 2.71) 

and  the  fluid  oscillator  equation  is 

w + w = (a[  - a'4)uc  wl  Da  - a'2w3/(ucD0)  + a\ucx\ID0  ( 2.72  ) 

where  ©*  is  the  vortex  shedding  frequency  given  by  co.,  = 2nSuc/D0  and  Strouhal 
number  S is  taken  as  equal  to  0.20. 

On  the  basis  of  the  results  of  various  vortex  experiments  for  subcritical  Reynolds  num- 
bers, Iwan  and  Blevins  choose 

a\  = a\la0,  o!2  - a2la0 , a\  = aja0 

a0  = 0.48,  a\  = 0.44,  a2  = 0.20,  a4  = 0.38  ( 2.73  ) 

Later,  Blevins  considered  the  spanwise  coupling  along  the  riser  because  portions  of 
the  riser  span  are  not  in  resonance  for  nonuniform  flow  (1990).  In  order  to  apply  the 
model  to  nonuniform  flow,  the  span  of  the  riser  is  divided  into  segments  that  are  within 
the  resonance  band  and  the  remainder,  which  is  outside  the  resonance  band. 
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Figure  2.5:  Control  Volume  Containing  a Cylinder  Shedding  Vortices 
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Define  a parameter  p (x3)  that  defines  the  spanwise  region  of  the  resonance  band: 


where  a and  P specify  the  lock-in  band  and  f„  = 2ti/co  , is  the  stationary  cylinder  vortex 
shedding  frequency.  Typically  a =0.6,  P = 1.4  for  lightly  damped,  large  amplitude 

motion.  The  parameter  p(x3)  will  step  between  0 and  1 along  the  span,  depending  on  the 
current  profile. 

The  evaluation  of  the  net  force  on  cylinder  under  the  consideration  of  the  spanwise 
coupling  is  obtained  from  the  modification  of  Eq.  (2.71)  and  is  given  by 

qi  = a4pwD0uc(w  - xi)p(x3)  - \pwucD0CvX\  (1  - p(x3))  (2.74) 

Since  the  model  derived  by  Blevin  and  Iwan  is  based  on  the  assumption  of  small 
amplitude,  this  model  is  applied  to  the  approximated  equation  of  motion  (2.61),  which 
gives 


From  above  equation,  we  can  recognize  that  the  riser  system  with  internal  flow  will 


1,  if  lock-in  band,  ( afs  </)  < P/,) 
0,  if  no  lock-in 


m,x i + 2mjV\x\  + mfV\x'(  + (Elxf)"  - wx,  - T ex'{ 

= a4  pw  Do  uc  ( w - xi  ) p(x3)  - \ pw  ucD0CdX\  (1  - p(x3)  ( 2.75  ) 


respond  resonantly  to  vortex  shedding  along  the  part  of  its  span  in  the  lock-in  band  and 
the  system  oscillation  will  be  damped  by  the  external  fluid  out  of  that  part. 


CHAPTER  3 

FINITE  ELEMENT  MODEL  OF  THE  SYSTEM 


This  chapter  formulates  the  numerical  model  of  governing  equations  and  boundary 
conditions  using  finite  element  method.  Variational  statement  for  the  boundary  value 
problem  is  introduced,  that  is,  the  weak  form  of  governing  equation  is  derived.  Sequen- 
tially, the  incremental  operator  is  applied  on  the  weak  form  in  order  to  derive  the  incre- 
mental equilibrium  equation  and  then,  Galerkin's  decomposition  method  is  performed  to 
construct  approximations  to  the  solutions  of  the  problem.  Finally  incorporation  of  the 
usual  cubic  shape  functions  yields  the  matrix  dynamic  equilibrium  equations  constructed 
in  terms  of  the  unknown  deformations  at  nodal  points.  In  addition,  the  finite  element  for- 
mulation for  vortex  induced  vibration  is  made  through  the  same  procedure  except  apply- 
ing incremental  operator. 

3.1  Variational  Statement  for  the  problem 
The  mathematical  model  of  the  present  problem  has  been  derived  in  chapter  2 and 
represented  in  vector  form  in  Eq.  (2.49).  To  produce  a variational  statement,  we  begin  by 
considering  Eq.  (2.49)  as  a two  point  boundary  value  problem 

m,  f +2w/Vir'  + #n/V?r"  + (Elr")" 

- [(Te  - ElK2)r/]/  - H(r'  x r")'  = q , 0 < j,  < Lx  ( 3.1  ) 

where 

A A A 

r = xi(ji)/  + x2(s\)j  + (s  + x3(si))k 

L\  - I (1+8/) cfc  , ( L = Riser  length  ) (3.2) 

JO 

This  is  partial  nonlinear  differential  equation,  with  respect  to  time  and  deformed  length, 
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of  fourth  order  with  variable  coefficients  defined  in  a domain  Q given  as  the  interval 
0 < ■?!  < L i . The  variable  coefficients  are  smooth  bounded  function  and  the  value  of  El  is 
always  positive.  Besides  Eq.  (2.49),  we  have  nonlinear  boundary  conditions  at  top  and 
bottom  of  riser,  i.e.,  sx  - 0 and  s\  = Lx.  Equations  (2.53)  - (2.55)  represents  the  top 
boundary  conditions  and  the  fixed-hinge  support  are  applied  at  the  lower  end,  that  is, 
those  boundary  conditions  can  be  given  as 

at  s i = 0 , r = 0 and  M = 0 

at  Si  = Zi , F]  = (TF)  t ■ i - K\  xiu 

F2  = (TF)  t j - K2x2u 

F3  = (TF)?  -k-K3x3u  + TTR  (3.3) 

and  the  torsional  moment  is  given  at  the  lower  or  upper  end. 

The  weak  form  of  the  mathematical  model  is  constructed  as  follows:  find  the  vec- 
tor function  r such  that  the  differential  equation  , together  with  the  boundary  conditions, 
are  satisfied  in  the  sense  of  weighted  average.  To  express  the  weighted  averages  of  dif- 
ferential equation,  we  first  construct  the  residual  error  function  vector  R of  Eq.  (3.1)  such 
as 

R(vi)  = mtr  + 2/m/Vi  r'  + m/V  \r"  + (EIr//)// 

- [(Te  - EIk2)!7]'  - H(r'  x r")'  - q , sx  e Q ( 3.4  ) 

AAA 

Multiply  R by  a sufficiently  smooth  test  function  vector,  r = xx  i + x2j  + x3  k,  defined 
over  the  entire  domain,  and  integrate  with  respect  to  the  deformed  arclength. 

J^Rrcfci  (3.5) 

The  weak  form  of  the  model  is  constructed  by  letting  the  above  integral  equal  to  zero. 
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The  finite  difference  form  of  the  model  is  given  by 

Nr  N r/j 


V f V P1 

2j  I R rds\  = 2j  I R rds\  - 0 

7=1  *^1  |=1  JO 


(3.6) 


where 


l 


'h 


lu  = I (1  + Bt)ds  , /,•  = the  length  of  i-th  element 


Substitute  Eq.  (3.4)  into  (3.6)  and  separate  the  integrand  for  each  subdomain  to  give 

Ch 

I [m,  r + 2/w/Vir'  + rrifV\x"  + (Elx")" 

JO 

- {(Te  - ElK2)r,}/  - H(r/  x r")'  - q ] • r ds\  - 0 (3  7) 


Integrate  by  parts  to  reduce  the  order  of  differentiation  to  introduce 

J^(EIr")"-r«*i  = (EIr,/)'.F|Jf-(EIr")-F/|Q1,+  ^ {Elx")  x" ds, 

Ch  if/] 

}0mfV*r"'?dsi  = mfvy  ■ r | q1'  - Jo  rrifNlr' -r' dsx 

Jq  ((Te  -ElK2)r'}'-F£*,  = {(Te  - EIk2)  r'}  • F | ^ ^ {(Te  - EIk2)  r'}  ■ f'dsl 


H(r'  x r")f -rdsi  = H(r'  x r")r  | - f H(r'  x x")x'dsi 


(3.8) 


Combining  Eq.  (3.8)  into  (3.7),  we  get 


Ch  _ Ch  _ pi  r/] 

JQ  mtx- rdsi  + Jq  2/W/Vi  x'  • r dsx  - JQ  mfN\x'  ■ x' dsx  + (El r")  • x"  dsx 

+ f {(Te  -EIk2)!*7}  • x'  ds\  + f H(r/xr//)F/^1=-/77/V?r/F|/n1  + (EIr//)?/ln 

Jo  Jo  0 ’0 

+ [-(EIr//)/+{(Te  - ElK2)r'}+  H(r'  x r")]  • F | g1'  + J q F cb\ 


(3.9) 
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From  Eq.  (2.14)  and  the  fact  stated  in  the  first  paragraph  of  Section  2.5.1,  we  have 

F =-(EIr//)/+{(Te  - ElK2)r'}+  H(r'  x r")  (3.10) 

and  also  we  have  the  following  relationship 

t x M = It  x (EIk£  + H?)  = -EIk/i  = -Elr"  (311) 


f*i 


Rewrite  Eq.  (3.9)  using  (3.10)  and  (3.1 1)  to  give 

[h  _ [h  _ r/i  _ n 

JQ  m,r  ■ rdsi  + 2/w/VIr/- rdsx  - mfV]r'  ■ r' dsx  + 

+ {(Te  -ElK2)r/}  • r' ds\  + J0H(r'  x r ")  ■ r'dsx  = {-/w/Vj r'  + F)  F | ^ 


(El  r")  r"ds\ 


(t  x M ) 


-ij.f 


q • r ds\ 


(3.12) 


Equation  (3.12)  represents  the  weak  formulation  of  the  boundary  value  problem  for 
each  subdomain  that  transforms  the  differential  equation  over  the  subdomain  and  bound- 
ary condition  at  the  both  end  into  a single  equation  in  which  all  of  the  features  of  the 
solution  and  the  discontinuous  data  are  intrinsically  present.  Moreover,  from  the  substitu- 
tion of  Eq.  (3.12)  into  (3.6),  we  can  define  the  same  variational  statement  over  the  entire 
system.  The  difference  is  that  the  boundary  conditions  represent  ones  at  the  both  end  of 
whole  domain  due  to  the  cancellation  of  the  boundary  conditions  of  each  subdomain. 

3.2  Incremental  Weak  Formulation 

In  previous  section,  the  weak  form  of  governing  equation  has  been  derived,  i.e., 
multiplying  both  sides  of  equation  by  a virtual  displacement  or  a test  function  and  inte- 
grating over  the  length  of  a finite  element.  And  then,  finite  approximations  on  the  weak 
form  yields  the  matrix  equilibrium  equation.  This  is  a usual  procedure  of  the  finite  ele- 
ment approximation.  However,  the  application  of  this  procedure  may  show  inaccuracy  in 


41 


solution  or  inefficiency  in  convergence  resulted  from  so  many  iteration  for  the  nonlinear 
terms  in  the  equation.  And  therefore,  the  application  of  incremental  operator  on  the  weak 
form  of  governing  equation  (Oden,  J.T.,  1972)  is  required  for  the  successful  iteration  and 
the  incremental  matrix  equilibrium  equation  can  be  obtained  from  the  incremental  weak 
form. 

The  incremental  operator  A has  the  meaning  of  a difference  operator,  i.e., 


where  x is  an  independent  variable  and  the  subscript  indicates  the  increment  number. 
The  operator  can  be  applied  to  a function / such  as: 


and  equation  (3.13)  can  be  rewritten  by  using  Taylor's  expansion.  For  example,  if  / is  a 
function  of  independent  variables  x andy , A/  can  be  expressed  in 


Now,  applying  incremental  operator  Aon  Eq.  (3.12)  with  independent  variable  r gives 


Ax  = x*+i  - xk 


A f - fk+  i - fk 


(3.13) 


A/  = ^Ax  + — Ay  + O (Ax2,  Ay2,  Ax  Ay) 


(3.14) 


mt  Ar  • r ds\  + 


Cl 


(EIAr  ")?"dsi  + 


+ JQ  {A(Te  -ElK2)r'  +(Te  -EIk2)  Ar'  } r'  dsi 


x r"  + r'  x Ar")  • r'dsi  = {-mfW i Ar'  + AF  } • r | lQl 


(3.15) 
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in  which 


A(Te  - EIk2  ) = ATe  - EIA(k2) 

Combine  Eqs.  (2.1),  (2.2),  and  (2.15)  into  (3.16),  we  obtain 
ATe  - EIA(k2)  = A(F  • r')  - EIA(  r"  • r")  = AF  • r'  + F • Ar'  - 2EIAr"  • r 
Substitute  the  result  of  (3.17)  into  (3.15),  we  finally  get 

Ch  _ [h 

J mtto-rdsl  + J 


(3.16) 


(3.17) 


, _ pi 

2ntfYi  Ar'  • rds\  - I nt/Vj  Ar'  ■ r'dsy 
0 JQ 


f/i  [h  f/i 

(EIAr ")  r"dsx  + Jq  (Te  -EIk2)  Ar'  • r'dsi  + 


(F-  Ar)(r/  • r')ds  i 


- (2Elr"  ■ Ar")(r'  ■ r')dsi  + H(Ar/  x r/,)-r/<fci 

f/j  , 

+ J H(r'  x Ar")  • r'  ds\  = {-/w/Vj  Ar'  + AF  } • r | J + EIAr"  • r'  | / ] 


f/l  _ f/i 

+ J Aq  • r ds\  - J (AF  • r/)(r/  • r')ds\ 


(3.18) 


This  is  the  incremental  weak  form  of  the  governing  equation  and  the  quantities  in 
front  of  integrals  are  considered  constants  for  each  element  so  that  those  terms  may  be 
factored  out  of  the  integrals.  In  order  to  represent  Eq.  (3.18)  into  components,  the  differ- 
entiating terms  of  the  position  vector  or  of  the  virtual  displacement  vector  with  respect 
to  time  and  deformed  arc  length  should  be  considered  and  are  given  by 

A A A 

r = xi  i + X2 j + (X3  + s)k 

AAA 

r = Xi  / + X2 j + x3  k 

.A  . A A 

r'  = x\  i + x2j  + x3k 
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r"  =x'(}  + x"j  + x"k 

I A A A 

r'  = Xj  i + x2j  + d\k 

A A A 

ir'  = Xj  / + x2j  + dq  k 

- . A A A 

r"  = Xj  / + x2j  + d^k 
k2  = xf  + x22  + d* 

Ar7  = Ax  j i + A x2j  + d<,k 

A A A 

Ar7  = Axj  i + Ax2j  + d(,k 
Ar"  = A x'l  i + Ax2 j + dqk 

AAA 

Ar  = Axi  i + A x2j  + Ax3  k ( 3.19  ) 

where 

. _ „i/2 

d\  — c j 

dq  = — C i C2 

dq  — — cj1/2  c3 

d<\  — c j 
d5  = - cj1/2  c4 
d6  = - (c73/2c2c4  + c71/2c5  ) 

dq  — — (Cj  C3C4  + Cj1/2C6) 

Cl  = 1 - X,2  - Xj2  = (x3  + 1/(1  + 8,))2 

C2  = XjXj  + XjXj 
C3  = X^Xj7  + XjXj 

c4  - XjAxj  + x2Ax2 
c 5 = AxjXi  +x/1Ax/1  + AX2X2  + x2Ax2 
c6  = A x[x'{  + XjAxf  + Ax2x2  + x2Ax2 


(3.20) 
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3.3  Finite  Element  Approximations 

The  construction  of  a finite  element  approximation  of  the  problem,  such  as  given  in 
(3.18),  is  based  on  Galerkin's  decomposition  method.  Here,  we  replace  (3.18)  by  a finite 
dimensional  problem  of  the  following  form  to  find  the  approximate  solution  vector  rh 
such  that 


f l\  _ f/i  _ r/j 

Jq  /w,Arh  • r^ds\  + 2/w/Vi  Af^  • rhtfci  - J mfV\ Ar£  ■ r'hdsx 


+ 


Jo  (EIAl'")T>.  + J'1  CT.-EI^ArJ-ri*.  + f'  (F.  Ar')(r'  •?')&, 


- r(2EIr?.Ar'')(r'  •?£)<*,  + J*1H(Ar,h  x rf)  •?'<*, 


+ J0H(ri  x4rh)r;<*i  = {-m/VfAri  + AF}rh  |(h  + EIAitf  • r£  | 'lf 


f/l  f/l 

+ Jo  Aq  rh*.  - J (AF-r'Xr' 


(3.21) 


The  nodal  values  of  the  approximate  solution  are  unknown  functions  of  time  and  consist 
of  not  only  the  deflections  at  the  nodes,  but  also  the  rotations.  In  other  words,  after  con- 
structing interpolation  functions  on  a suitable  mesh,  we  take  the  approximate  solution  to 
be  of  the  form 


rh  = xih7  + x2hj  + (x3h  + s)k 


xjh(si,t)  = E xM (t)Ni(sl) 
/=  i 


and 


j = 1,2,3 


( 3.22  ) 
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where  x7i , xfl  and  xJ2,xj4  represent  the  deflections  and  the  rotations  at  node,  and  Nj  is  the 
basis  or  interpolation  function. 

As  an  interpolation  function,  Hermite  polynomials  are  constructed  by  matching  to- 
gether element  shape  functions  in  the  usual  way.  Let  ft  denote  a master  element,  with 
coordinate  - 1 < £ < 1 . Note  that  the  specification  of  xyh  and  x'jh  at  the  endpoints  of  an 
element  involves  four  conditions.  Since  a complete  cubic  in  £ contains  exactly  four  pa- 
rameters, we  can  construct  cubic  shape  functions  Af;  which  exhibit  the  necessary  proper- 
ties: 


Ni  ( - 1 ) = 1,  N i ( 1 ) = 0,  dNi(-  1 )/</£  = dNi(\)/dt,  = 0 
N3(-l)  = 0,  #3(1)=1,  dN3(-l)/d$  = dN3(l)/d$  = 0 
N2(-l)  = N2(l)  = 0,  dN2(-\)/dt  = 1,  dN2(l)/dt  = 0 
N4(-\)  = #4(1)  = 0,  dNA(-l)/dZ,  = 0,  dNA(l)/d$=l  (3.23) 

These  conditions  define  a set  of  unique  local  cubic  Hermite  polynomial  shape  functions, 
#i(t)  = J(2-3£+53),  &($)  = i(2  + 35-^3) 


= \(-i-^+e+e)  0.24) 

which  are  illustrated  in  Fig.  3.1. 

The  restriction  x^  of  a approximate  solution  xju  to  an  element  Q,  between  nodes  / 
and  / + 1 can  be  written  as 


xjh(*i,0  = *jh(5(*i),0 


( 3.25  ) 


where 


2si-(su+si+1) 

5(*l)  = " 


‘Hi+l  Xjj 


(3.26) 
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and  xjh  is  the  approximated  solution  referred  to  the  master  element, 

= txj,(t)N,a)  (3.27) 

/=  1 

The  element  shape  function  Nj  is  obtained  from  (3.24)  and  (3.26)  as 


Ni(Si)  = #,(€($i))  (3.28) 

Combining  corresponding  shape  functions  Nj  at  node  i determines  the  global  basis  func- 
tion NJ-  at  node  i (Fig.  3.1). 


4 


Figure  3.1:  (a)  Hermite  Cubic  Shape  Functions  on  a Master  Element  and  (b)  Correspond- 
ing Global  Basis  Functions  at  a Node  i in  the  Mesh. 
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Having  selected  the  basis  functions,  the  incremental  form  of  a approximate  solution 
and  the  test  function  can  be  written  similarly 

4 

Axjh(ji,/)  = ZA xji(t)Ni(si) 

7=1 

4 

**($1,0=  (3.29) 

7=1 

The  substitution  of  (3.29)  into  (3.21)  with  the  consideration  of  (3.19)  introduces 
the  series  form  of  the  incremental  weak  form.  Manipulating  the  series  form  and  writing 
term  by  term, 


r 

Jo 


fl  1 

m,Arh-rh<*i  = /W7(Axihxih  + Ax2hX2h  + AxshXsh)-*! 


= f lm,{  2 ( E ticjiNi)(  Z AxjiNi)}dsi 

Jo  j=  1 7=4  7=4 

3 4 4 r/j 

= Z [ Z { Z ( mtNiNjdsi)A3tu}xu] 
fc= i t=i  m Jo 


t 


2m/Vi  Arh  ■ rhcfti  = I 2m/Vi(Axlhxlh  + Ax2hx2h  + 

0 J0 

f 1 2mfVi  { Z ( Z tejiN'i  )(  Z XjiNi  ) + d6  Z x3 ,Nt  }ds] 

JO 


Ch  ± 4 

'0  >1  7=1 

2 4 4 f/] 

Z " " 

/c=l 


4 

Z 

/=1 


4 

Z 

7=1 


(3.30) 


2 4 4 f/l  4 r/j 

Z [ Z { Z ( 2mf\\NiN,jds\  )Ax^>xu]  + S ( hn/VidsNidsy  )x3i 

k=  I 7=1  >1  J0  Pi  JO 


(3.31) 


f ’1-»/V?Ar,h  rh^,=  Z [ Z { Z ( )Ax#}5c«] 

J0  A=1  7=1  7=1  JO 


+ 


4 Cl  1 

Z ( I -mf\\  dil/jdsx  )x3i 

7=1  J0 


(3.32) 
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f*  El  Aif  = t [ t { Z ( f^EI  Nf'N''dsl  )A xkj}xki] 

J0  *=1  i=l  >1  Jo  y 

4 r/j 

+ Z(JQ  EI  d.-jN'l ds\  )x3/  (3.33) 


f ^(Te  -ElK2)Ar'  •«•{>!  = Z [ £ { £ ( fVe  - El k2  ) A^ A^,  ) Axy  } xw  ] 

JO  *=1  r=l  >1  JO 

4 f/j 

+ Z(  (Te  - EIk2)  dsN/idsi  )x3;  (3.34) 

<=  l JO 


7l(F- Ar£)(r£  ■!■£)<&! 
0 


A*lh  F2AX2h  + F3£Aj  ) (*lh*lh  +X2h  X2h  + d\  Xj^  ) ds 


4 4 r/j  r/j 

= Z Z {(  J^ix'^A^OAxy  + (JoF2x/lh^A5*,)Ax27}xu 

44  r/j  r/| 

+ Z Z {(  Jo  Flx'2hN/iN/jds1)Axlj  + ( Jo  F2x'hA^7V'^1)Ax27  }x2/ 

4 4 f/i  f/i 

+ Z Z {(  Jo  F ! d\  A^c&OAxy  + ( Jo  F 2di  N/iN/jdsl)Ax2j}x3i 

4 f h X f h 

+ Z ( F3  J5Xih  <*1  )xi;  + Z ( F3rf5X2h  JV|fl!si)jc2; 

*=i  JO  *=1  JO 

4 r/] 

+ Z ( F 3 d$  d\  A/y  ds  1 )x3/ 

j=i  Jo 


(3.35) 


f/i 

Jo-2EI(r'/  Ar'/)(r' •?')*, 
fZi 

- -2EI(xlhAxIh  + x2hAx2h  + d^di)  (xlhxlh 


+ x2hx2h  + ^i*3h)<*i 
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- 5 | {(  {0'-2El4<hA^'*,)Ax,,  + ( ds^tXiAxu 

E E {(  J1-2EIx'ah4A;^*i)A*u  + ( j ]-2EI x'2hx"hN/iN',dsl)Ax2j}x2i 


+ 


+ 


/=  i >i  Jo 
f/l 


E E {(  f1  _2 EI xfh diN/jN,j' ds\)Axy  + ( f^EIx"  * JV^'ds,  )Ax2,}x3, 

j=i  ^=i  •'U  •'U 

4 r/i  4 r/! 

E ( -2EIfli,3rf7*,|h^«fci)xi,  + E ( -2EId3 d7x'2bN'ids i)x2i 

j=l  JO  f=i  Jo 


f=i  JO 
f/l 


4 f/l 

+ E ( I -2EIJ]  d3dq  iJjdsx  )x3/- 
#-i  Jo 

f/l 

J H(Arh  x rh  ) ‘ rh^l 


(3.36) 


ri 

— I H{(£/3Ax2h  — £/5X2j1)xlh  + (i/5Xlh  - (af3Axlh)x2h  + (x2hAxlh  — xlhAx2h)x3h  }£&i 

= E E(  f1Hrf3iVjiy;*1)Ax^li+  E E ( j*^1  -Ha/3  N,jN/jdsi)  Axi/x2/ 

r=l  >1  JO  i=  1 >1  JO 

4 4 f/l  44  f/l 

+ E E ( Jo  Hx2h N,iN,jds\) AxijXn  + E E ( Jq  -Hx"  ) Ax2Jx3, 

+ E ( \l]-Bd5xiN'idsl)xu+  E ( pHrfjxJ'h^rfsOxa  (3.37) 

i=i  JO  i=i  JO 

J0  H(rh  x Ar b)  r'hdsi 
Cl  i 

— H{(6?7X2h  — <7l  Ax2h  )xlh  + (<7iAxlh  — d-jX^^X^  "1"  (XjhAx2h  — ^2hA^ih  )x3h  }ds\ 

4 4 f/l  44  f/l 

= E E ( -Hflf,  N/iN''dsl)AxyXu  + E E ( N* N1' dsx  ) Axyx2/ 

/=i  >i  JO  /=i  >i  Jo 

4 4 f/l  44  (*/] 

+ E E ( I Hx'lhN,iN''dsl)Ax2jx3i+  E E ( -Hx2h  N/iN,j/  ds\  ) Axi7x3; 

j=i  >=i  •'O  r=i  >=i  •'O 
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+ 


4 cl\  4 r/j 

I (JoHrf7x5h^*i)x1,+  E (Jo-Hrf7x'hJV'£&1)jc2/ 


M JO 


f/l  4 f/l 

Aq  ■ r h i = (Aqixih  + Aq2x2h  + Aq3x3h)^i 

Jo  Jo 

3 4 f/i 

= 2 2(1. 

1 t=1  0 


Aq*V,c&i  )xki 


f/l 

— — 1 (AFixlh  + AF2X2h  + AF3£/i  ) (xlhxlh  + *2h*2h  "*"^1*311)^ 

2 4 f/i 

= E [ £ { J -(AFix'h  + AF2x2h  + AF3(/i  ) x'faN'idsi  }xfe] 

/o=l  f=l  * vJ 

f/l 


+ 


4 f/j 

2 { I -(AFix^b  + AF2X2h  + AF3Ji  ) d\  N'j  ds\  }x 
/=i  JO 


(3.38) 


(3.39) 


( 3.40 ) 


( - m-,  Vi  Ar£  + AF  ) • r h | lQli  = - /w,-  Vj  Ar£  • rh  + AF  • rh|  ^ 

= -ntiV  1 (Ax'lhxih  + Ax2hx2h  + d5x3h)  +(AFixIh  + AF2x2h  + AF3x3h)  | ^ 

= [ -mi V\  {E  ( E E N'j Nj AxiaXkj)  + d5(t  Mx3;) } + E { AF*(  tv,-  x«  ) }]  1 7» 

*=1  >1  i=  1 A=1  r=l  1 0 

= -mtV\  {(Axi4Xi3  - AX12X11)  + (Ax24^23  - AX22X21)  + (d5(lu)x 33  — d$ (0) X 3 1 ) } 

3 

+ E { AFjfc(/i,-)x*3  - AF*(0)xti  } ( 3.41  ) 

k=  1 


EI4r".?'h  |'"  = H(A41*'h  + A4.i'h  + rf7x'„  |'» 

= [ HI  { i(t  t/V;'/V'Ax1,x*))  + d7(t/V'S3,)}l|(1' 

k=  1 /=!  /=!  f=l  ” 


= (AMi(/i,)xi4  -AMi(0)xi2)  +(AM2(/i,)x24  -AM2(0)x22) 
+ ( AM3(/ u)x34  -AM3(0)x32) 


( 3.42  ) 
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3.4  Matrix  Dynamic  Equilibrium  Equation 
In  previous  section,  The  finite  element  representations  for  each  individual  compo- 
nent of  the  incremental  weak  form  have  been  given  in  the  series  form  of  Equations  (3.30) 
- (3.42)  through  the  incorporation  of  Hermite  polynomials  selected  as  the  shape  or  basis 
functions.  A total  of  12  degrees  of  freedom  per  element  (six  for  each  node)  have  been 
used  for  the  representation  of  the  deformation,  namely  the  deflections  and  their  deriv- 
atives (rotations). 

To  determine  the  deflections  and  the  rotations  at  node  (x7i , xj3  and  x# ,x7-4,  j =1,2,3) 
of  coefficients  that  will  characterize  the  approximate  solution  in  (3.27),  we  substitute 
(3.30)  through  (3.42)  into  (3.21)  to  obtain  the  condition 


mtNj Njds i ) A % }xki] 


+ 


2mfViNiN,jdsi  )Axkj}xki] 


+ 


2mf\\d6Njds\  )x3j- 


4 f^i 

-mf\1lN,iN,jdsx  )Ax#}xh]  + £ ( I -mf\\  d^Nlids\  )x3l- 

f=l  “O 

Elrf' rfj' dsx)bxkj}xki]  + £(  f \\d1N,i,dsx  )x3, 

i=  1 »0 

(Te  - ElK2)A^7V'iA:1)Ax^}xfa] 


4 Ch 

+ £(  (Te-ElK  2)d5N'idsl)x3i 
i=  1 JO 

4 4 f/ 1 

+ £ I{(  + ( I ¥2x[hN/i  N/jdsi)Ax2j}xu 

r=l  >=1  •'0  •'O 

4 4 r/ 1 r/i 

+ I I {(  J0  F lx'2hN'iN,jdsl)Axlj  + ( Jo  F 2x,2hN/iN,jdsl)Ax2j}x2i 
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+ 


+ 


4 4 f/]  f/j 

E E {(  Jo  F i d\  N'iN/jds\)Axij  + ( Jq  F2dx  N/iN/jdsx ) Ax2j  }x3; 

4 f/i  4 f/i 

5 ( I F3 d$ TV'  *i)x1(-  + E ( I F3*x'2h  N',  dsx)x2i 

j=  1 JO  r=l  JO 

E ( F3 a?5 d\  TV'  ds i)x3i 

i- 1 v ( ) 


(=1  JO 


+ 


+ 


4 4 r/ 1 r/] 

E E {(  Jo  -2  EI  x jhx  i'h  N'jbfJ  ds\)^x\j  + ( Jo  -2EIx/lhx'/h^A5/cfe1)Ax27}xH 

44  r/ 1 I*/ 1 

E E {(  JQ  -2EIx/2hx'/hTV'TV''^i)Axu-  + ( JQ  -2EIx'hx'/hTV'A^/cfe1 ) Ax2y- }x2/ 

E E {(  f -2EI  x'{hd\  N'jN'j  dsx  )Axi7-  + ( I*  -2EIx2h^i N,iN/'dsx  ) Ax^-  }x3j- 

r=l  j=  1 JO  JO 

4 r/i  4 r/ 1 

E ( I -2EIc/3  J7XihTV'fl5s'i)xi/+  E ( I -2EIaf3iT7X2hA^cTsi)x2; 

f=i  JO  <=  i JO 


+ 


+ 


+ 


+ 


+ 


+ 


4 f h 

E ( I -2EIcTi  c/3  £^7  A^c&i  )x3;- 

/=i  JO 

4 4 r/j  4 4 r/] 

E E(  I HdjN'jN'jdsx  ) Ax2;Xi;  + E E(  I -Hd3  l/jN'jdsi  ) AxyX2; 

/=i  >1  Jo  m >1  JO 

4 4 r/j  44  r/i 

E E ( Hx2hN/iN,:dsi)AxljX3i+  E E(  -Hx'{hN,iN,Jdsx  ) Ax27x3/ 

;=1  >i  Jo  »=i  >1  JO 

4 f“T  1 4 | ■/] 

E ( L -Hrf5X2h^(*i)xii  + E ( H<V5x^hA^i*i)x2i 

F=1  JO  j=  l JU 

E E ( L N'jN'j' dsi)Ax2jXU  + E E ( [ Hfi?i  A^jVfc&i ) Axyx2/- 

r=l  >1  J0  M >1  J0 

4 4 r/i  4 4 r/] 

E E ( Hx'lh  N'iN1' dsx)^x2jX3i  + E E ( Jo  -Hx'h  A^A^i)Axi7-x3/ 

4 f/l  4 f/l 

E ( HtT7  x2h  ds  i ) x if  + El  ( -Ud7x\hN/idsx)x2i 
7=1  JO  7=1  JO 
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- Z Z ( AqkNidsi  )xki 

A=1  f=  1 «0 

2 4 Cl\ 

+ I [ I { I -(AFjx'h  + AF2x'h  + AF3rf,  )x'khN'idsl  }xki ] 

A=1  f=l  •'U 

4 r/i 

+ Z { -( AFix^h  + AF 2*2ii  + AF3^i  )dxN/idsi  }x'v 

j=i  *’0 


-mfV]  {(Axi4Xi3  - AX12X11)  + (Ax24*23  - Ax22X2i ) + (d$(l i;)  x33  -d5(0)x3i)  } 


3 

+ Z { AF*(/i/)x*3  - AF*(0)x*i  } 

A=1 


+ (AMi(/i,)x,4  -AMi(0)xi2)  +(AM2(/i/)x24  -AM2(0)x22) 


+ (AM3(/i;)X34  — AM3  (0)  X 32  ) 


( 3.43  ) 


Factorizing  the  coefficients  of  test  function  (x«,  xw,  and  x*3,  k = 1, 2, 3 from  Eq. 
(3.43),  we  may  have  the  following  compact  form 

3 4 _ 4 3 4 

Z [ Z xkj  { z (Mij  Ax#  + Cij  Ax#  + Ky  Axkj ) } ] = Z Z(  A/jy  + AFjy)x;t  (3.44) 

A=1  »=i  >1  lc=  1 /=l 

Since  the  x«  are  arbitrary,  (3.44)  represents  12  independent  equations  to  be  satis- 
fied by  the  Ax#  rather  than  the  single  equation  it  may  appear  to  be.  In  other  words,  the 
coefficients  of  test  function  may  be  taken  out  by  the  natural  choices  for  the  set  of  xki. 

For  an  instance,  for  the  choice  of  x,,  = 1,  x12  = x13  = xM  = = x33  =x34  = 0 we 

have 

4 4 

Z (MijAxy  + C\j Ax#  + K\j Ax#) }]  = Z (A/),-  + AFi;)xi;  (3.45) 

>1  t=i 

Continuing  this  way,  we  arrive  at  the  system  of  12  linear  independent  equations  in  the  12 
unknown  coefficients  of  approximate  solution,  namely  A xkj,  k = 1,2,3,  and  j = 1,2, 3, 4. 
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Writing  it  in  matrix  form, 


[My] 

[0] 

[0] 

{A*  17} 

[Cij\  [0]  [0] 

{Mi,} 

{ 0 } 

[0] 

[Mv] 

[0] 

{Ax2j} 

► + 

[0]  [Cy]  [0] 

(Ax  27} 

• + • 

{0} 

[0] 

[0] 

[My] 

{Axy} 

[0]  [0]  [0] 

{Ax3y} 

ify} 

{fk4j  + fk5j  +fk6j  + fklj\ 

Wj  +f/j} 

[Cey]  [0]  [0] 

{A*  17} 

+ • 

\fk4j  + fk5j  + fk6j  +fk7j} 

* — < 

► + 

[0]  [Cey]  [0] 

(Ax  27} 

\fklj  +fk2j  +/k3j  + fk4j  + fkSj\ 

. Wj+ffj\ . 

[0]  [0]  [0] 

{AX37} 

{0} 

{AF,;} 

+ • 

{0} 

• + • 

(AF2j}  ■ 

{/cej} 

{AFjj} 

(3.46) 


where, 

My  = 1 

f/[ 

| mtNjNjdsi  , 

h = f 

h 

2mf\\deNjds\ , 
3 

fk\j  = j 

|*i 

|Q 

fhj  = j 

f/l 

ElrfvM'flfai, 
'0  J 

= j 

f 1 (Te  - EIk2  ) d<,Nijds\ 

Kmj  = 

[/if2x/Uia'a;^1, 

Cy  = ^l2mfWiNiN,jdsx 

[h 

Kvj  = J o-mfViN'iN'jdsl 

[h 

K2ij  = JQ  El  N't  rfj'dsx 

fll 

Kjjj  = (Te  -ElK2)N'iN'jdsx 

Jo 

cl  i 

K"  = JQ  F, x[hN/jN/jds\ 

K\)j  = J]F  lx'2hN/iN'dsi 


55 


r^22 

A4  ij  - 

^ = 1 

h 

Fi  d\N,jN,jds\ 
0 J 

zX32 

A4  ij  ” 

A = | 

) F 3 d$  N'  dsi 

fm  = \ 

f/i 

|0  F3^5X2h  TVj  ds  1, 

A = J, 

h 

F3  ds  d\  Nj  ds\ 

II 

J’o'-2EIx',h<^<*1, 

= J 

*1 

0 -2EIx;„x"  f/^ds, 

II 

(*/l 

*3  ■ j 

•h 

o-2EIxi,x"  l/rfdsi 

II 

— 

m «o 

f/ 1 

Jo-2Ei*Jlrf1A;jyf<&1, 

*3  - J 

0-2EI  x"  d.W'w''*, 

fkSj  - J 

r/i 

1 -2EIc/3C/7X/lhM<&i, 
0 

A = j 

'll 

-2Elds  d-j  x'2h  Njds\ 

fk5j  = j 

-2  El  d3  di  l/:  ds\ , 
0 J 

r^l2 

A6//  “ 

f/j 

Jq  Hd3N'lN'Jdsl 

t^21  _ 

^6  ij  ~ 

f '-HrfjMiVlffei, 
'0 

*3  = J 

f/i 

|Q  H x'i^N'ds, 

*3  = J 

f/l 

/k6/  = 

f/i 

| ~ H ds  x2h  rfj  d.s  i 

A = j 

'l\ 

Hd5x"hN'dsu 

0 

1^12 

A7t/  “ 

f/i 

Jo  -H dx  N'iN''dsl 

^ = j 

f/i 

UdlN,iN//dsu 

|Q 

rX31 

A7t/  “ 

f/i 

]0Hx'lhN'iN',ds1 

*5  “ J 

f1  ~Hx2h  l/jf/j  dsi, 
in 

/k7/  ~ 

j*  Hd7X2hN/jds\ 

Jo 
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qj  = I Aqi  Njdsi 
0 


fi,  = y^-HdW^dsu  q)  = f' 

r/|  f/i 

q)  = Jo  &(\iNjds\  q]  = Jo  Aq \Njds\ 

[1\ 

f/j  = JQ  “(AFi4  + AF2x'h  + AF3*  ) x'lh N'ds, 

fij  = ^-(AFjx,,  + AF2x'h  + AF3^i  ) x'2hN'Jdsl 
•'0 

rl\ 

% = -(AFiXjh  + AF2x'h  + AF )d]N/jds1 

JO 


{AFi,}  = 


-AF  i (0) 
-AM , (0) 
AF,(/n) 
AM ,(/]/) 


{^2 ;}  = 


-AF2(0) 

-AM2(0) 

af2(/„) 

am2(/„) 


{My}  = 


( 3.47  ) 


-AF3(0) 

-AM3(0) 

af3(/„) 

am3(/i;) 


[Ceij]  = -JW/V? 


0-100 
0 0 0 0 

0 0 0 1 

0 0 0 0 


{ feej } - W/Vj 


-ds(0) 

0 

-d,(lu) 

0 


(3.48) 


In  equation  (3.46),  the  subscript  i and  j run  from  1 to  4 and  a total  12  degrees  of 
freedom  per  element  are  used  as  shown  in  Figure  3.2. 


X* 

-> 


Figure  3.2:  An  Element  of  Riser  Showing  Forces  and  Displacements  at  the  Node 


57 


Also  we  may  represent  equation  (3.46)  in  the  most  compact  form: 

[M]{A*}  + [C]{Ax  } + [AT]{Ax  } = {A/}  + {AF  } (3.49) 

where, 


[Mij]  [0]  [0] 

[Cij]  [0]  [0] 

[M]  = 

[0]  [M}j\  [0] 

[< c ] = 

[0]  [C,j\  [0] 

[0]  [0]  [Mjj\ 

[0]  [0]  [0]  J 

"[* 

ij  +K2  ij  +KVj  +K\)j  +K\)j ] 

\_K\l  +K\fj  +KX1  +Kl7fj 

[0] " 

[K]  = 

21  . 1^21  , l/l\  yl\ 

_4ij  +A5  ij  + A6  ij  +Klij_ 

[_Kuj  +K2jj  + Kjij  +K%j  +K%\  [ 0 ] 

1^31  y31  y3 1 y3 1 

_ 4i/'  + A5  ij  +A6(7  +A7(7_ 

[K%  +Kll  +Kll  +K% 

[°] 

[Cejf]  [0]  [0] 
[0]  [CeV\  [0] 
[0]  [0]  [0] 


Wi  +f/j  1 

[fk4j  +fk5j  + fk6j  + fklj) 

f > 

{0} 

' {0} 

(A/}  = - 

> _ < 

\fk4j  +fk5j  +fk6j  + /k7/} 

> < 

{0} 

■ + • 

{0} 

iW+jg}  J 

\fk\j  + fklj  +fk3j  + fk4j  +fk5j] 

{fcj) 

{ fceJ } 

{AF} 


(AF  V> 

(AFV)  • 
{**v} 


(3.50) 


In  (3.49),  the  matrix  [ M ] is  a consistent  mass  matrix  which  is  identical  to  the  one 
derived  for  small  deformation  models  and  includes  the  added  mass  of  the  riser  and  the 
mass  of  the  drilling  fluid.  [ C ] is  an  anti-symmetric  matrix  due  to  the  coriolis  force  re- 
sulting from  the  internal  flow.  Furthermore  [ K ] is  the  time  and  deformation  dependent 
stiffness  matrix  and  consists  of  8 components:  [ Kx  ] is  due  to  the  centrifugal  force  from 
the  internal  flowing  fluid,  [ K2  ] is  the  bending  stiffness  matrix,  and  [ AT3  ] is  the  geomet- 
ric stiffness  matrix  which  implicates  the  nonlinear  coupling  between  axial  tension  and 
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bending.  The  previous  three  matrixes  are  symmetric,  while  [ KA  ],  [ K5  ],  [ K6  ] and  [ Kn  ] 
are  non-symmetric,  full  matrixes.  [ K4  ],  [ K,  ] are  due  to  the  variable  longitudinal  ten- 
sion and  [ AT6  ],  [ AT7  ] are  due  to  the  torque  H.  Finally,  [ Ce  ] is  the  centrifugal  force  com- 
ponent at  the  both  end  of  element.  As  forcing  components,  { / } is  the  deformation 
dependent  equivalent  nodal  force  vector  and  { F } is  the  internal  force  and  moment 
vector. 

As  stated  in  section  2.5.1,  utilization  of  the  kinematic  constraint  on  unit  tangent 
vector  leads  to  the  reduction  of  the  degrees  of  freedom  per  element  from  12  to  8,  that  is, 
the  vertical  degrees  of  freedom  can  be  computed  in  terms  of  the  lateral  ones.  Thus,  we 
may  need  to  evaluate  the  lateral  degrees  from  Eq.  (3.46).  Manipulating  (3.46)  through 
the  partition  of  matrix  and  rewriting  only  in  terms  of  lateral  degrees, 


[Mij]  [0] 

[0]  [M0] 


{My}  1 I"  [Cv]  [0] 
{Axy}  \ |_  [ o ] [Cij\ 


{Axy}  | 

{Ax2j}  J 


+ 


[K\ij  + Kjij  + Kyj  + K4]j  + K\)j  - Cej\  \_K\l  + + Kly  + 

[K\)j  +K%  +K%  +K%\  [KUj+K2iJ+K3ij +K%+K]l  - CeiJ] 


W +f/j 1 

M) 


\fk4j  + fk5j  +fhj  +fwij) 
\fk4j  + fkSj  +/k6j  +/k7/} 


> + 


{AFb} 


{AXy}  1 

{Axy}  J 


(3.51) 


For  the  vertical  degrees  of  freedom,  the  matrix  equilibrium  can  be  written  as: 

im{tey)  + fa)^[K%+K%+K%+K%~\  [k^j+k^j+kSj+k^  ]|  j 


+ {fk\j  + fk2j  + /k3/  +fk4j  + fkSj]  - Wj  +f/j  ) + {/cej  > + (AF37  } 


(3.52) 
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Since  the  vibrations  in  vertical  direction  are  not  important,  we  can  ignore  the  iner- 
tia term  in  (3.52)  and  obtain  the  internal  forces  and  moment  from  (3.52)  as  following: 

{AF37  } = - {/cej  } + {/kly  +fuj  +f^j  +/k34y.  +/ky  - \q]  +fyj  } + {fcj} 

+ [ [K\]j  + K\]j  + Kl)j + ^7,7 ] +K%  +k%  +KTij]  ]|  1^1  J (3.53) 

Having  calculated  the  matrixes  and  equations  describing  our  approximation  over 
each  finite  element,  the  next  step  is  to  assemble  the  equations  describing  the  approxima- 
tion on  the  entire  mesh  by  adding  up  the  contributions  to  these  equations  furnished  by 
each  element.  This  assembling  procedure  is  omitted  because  it  is  well  known  and  refer- 
enced in  several  books  and  papers. 

3.5  Finite  Element  Model  of  Linear  System 

This  section  formulates  the  numerical  model  of  linear  governing  equations  derived 
under  the  assumption  of  small  deformation,  that  is,  transverse  vibration  and  in-line  dy- 
namic response  of  current-vortex  model  are  modeled  using  finite  element  method.  The 
same  procedure  as  discussed  in  previous  sections  ( section  3.1  through  3.4  ) is  applied 
except  the  utilization  of  the  incremental  operator.  The  application  of  the  incremental  op- 
erator is  not  necessary  in  this  problem  because  the  inaccuracy  of  solution  or  the  ineffi- 
ciency in  convergence  is  not  expected. 

The  mathematical  model  of  current  induced  vibration  can  be  obtained  simply  by 
combining  Eq.  (2.69),  (2.72)  and  (2.75): 

mtx  1 + 2m f\ \x\  + mfV \x'[  + (ELcf)"  - wx[  - Tex'{ 

= a4pwD0uc(w  - x{)p(x3)  - l-pwucD0CDxi(\  - p(x3)  (3.54) 

w + = (a[  - a'4)uc  wl  D0  - a'2w3/ (ucD0)  + a'4ucxx/D0  ( 3.55  ) 

mtx2  + 2m/Vix'2  + mjV\x2  + (JBAx2)''  - wx2  - = q2 


(3.56) 


60 


As  we  discussed  in  chapter  2,  those  three  equations  represent  the  transverse  vortex 
induced  vibration,  the  fluid  oscillator  equation  about  the  hidden  flow  variable  w,  and  the 
riser  vibration  in  the  current  direction.  The  forcing  term  q2  in  the  right  side  of  Eq.  (3.56) 
can  be  derived  by  combining  Eq.  (2.38),  (2.39),  and  (2.48),  and  only  taking  the  drag  term 
due  to  current  and  given  as, 


The  weak  form  of  the  equations  are  produced  in  multiplying  the  residual  by  a suffi- 
ciently smooth  test  function,  integrating  the  product  by  parts,  and  equating  the  result  to 
zero.  The  details  in  derivation  are  given  in  Appendix  B and  only  the  results  are  written 
here  as  followings: 


q2  = - jP>vCdD0  |x2  - uc|(x2  - uc) 


(3.57) 
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(3.60) 


Implementation  of  Hermite  polynomials  as  basis  functions  yields  the  matrix  dy- 
namic equilibrium  equations  constructed  in  terms  of  the  unknown  deformations  at  node. 
For  the  construction  of  the  matrix  equation,  the  same  processes  applied  in  section  3.3  and 
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3.4  are  used.  The  details  in  derivation  are  given  in  Appendix  B and  the  resulting  matrix 
equations  are: 

[Mij ] + [Cij  + CfiJ ] (x,//)}  +[Tj-  Ceij  + Kij ] (x,//)}  = [Fij ] {wjit)}  ( 3.61  ) 

Wij ] {%(0}  + [Cij ] + [Kjj ] {wj{t)}  = [Fij ]{Xlj}-{fj}  ( 3.62  ) 

[Mij  ] {*2/0}  + [Cij  ] (X2/0)  + [fij  - Ceij  + kij  ] {x2/0}  = {fj  } ( 3.63  ) 
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Eqs.  (3.61)  and  (3.62)  represent  the  vortex  induced  vibration  which  characterizes 
self-excited  oscillator.  Once  the  forcing  due  to  well-formed  vorticities  acts  on  riser,  the 
velocity  of  riser  itself  act  again  to  fluid  particle  as  a forcing.  Eq.  (3.63)  presents  in-line 
vibration  due  to  current.  The  characteristics  of  the  coefficient  matrix  in  (3.64)  will  be 
discussed  in  section  6.3. 


CHAPTER  4 

DEVELOPMENT  OF  A FINITE  ELEMENT  PROGRAM 


In  this  chapter,  we  discuss  the  development  of  a Finite  element  code  for  the  solution 
of  our  two-point  boundary  value  problem.  Some  considerations  of  code  design  are  pres- 
ented and  overall  structure  of  the  code  is  outlined.  However,  it  should  be  noted  that  the 
time  step  integration  and  the  predictor-corrector  scheme  for  the  convergence  in  solution 
due  to  the  nonlinear  terms  (subroutine  NEWMRK)  is  not  described  in  this  chapter  but 
treated  in  chapter  6. 

4.1  Description  of  CODE 

We  now  describe  the  features  of  our  finite  element  code,  called  CODE  which  is 
designed  to  perform  the  finite  element  analysis  of  our  two-point  boundary  value  problem. 
Details  on  the  construction  of  CODE  and  the  various  subroutines  it  contains  are  dis- 
cussed in  the  remainder  of  this  chapter.  As  viewed  in  Figure  4.1,  CODE  consists  of  the 
two  functional  units.  Since  the  finite  element  modeling  is  defined  by  the  input  data,  the 
reading  and  generation  of  the  arrays  of  data  are  performed  by  the  preprocessing  unit  of 
the  program  which  consists  of  several  subroutines.  The  actual  finite  element  calculations 


PREPROCESSOR 

PROCESSOR 

Read  input  data. 

Set  up  and  solve 

Define  parameters 

equations  for  nodal 

and  data  arrays 

point  values  of 

solution. 

Figure  4.1:  Functional  Units  of  the  Finite  Element  Code 
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are  carried  out  in  the  processor.  In  this  part  of  program,  element  matrices  and  vectors  are 
calculated  and  assembled;  boundary  conditions  are  imposed,  and  a set  of  global  equations 
is  solved  for  the  nodal-point  values  of  solution. 

The  flowchart  in  Figure  4.2  shows  each  subroutine  in  CODE  and  its  relation  to  oth- 
er routines.  Routine  names  are  given  in  boxes  and  a brief  description  of  their  function  is 
given  in  Figure  4.3.  The  data  in  CODE  are  organized  into  several  labeled  COMMON 
blocks,  whose  contents  are  summarized  in  Figure  4.4.  As  a general  rule,  the  data  are  de- 
fined in  routines  called  by  the  preprocessor  PREP  and  used  through  CODE. 

The  main  program  serves  only  to  call  the  routine  SETINT  that  initializes  numerical 
integration  data  (see  Section  4.2)  and  the  two  main  units  of  CODE  as  shown  in  Fig.  4.2. 

4.2  Element  Calculations 

In  the  center  of  any  finite  element  program,  the  calculation  and  assembly  of  the 
element  matrixes  and  vectors  are  performed.  These  calculations  consist  of  the  integration 
indicated  in  the  variational  statement  and  are  carried  out  in  subroutine  ELEMM  and 
ELEMV  in  conjunction  with  SETINT,  SHAPE,  GETCON,  and  GETMAT. 

In  order  to  compute  the  quantities  of  the  integration  in  (3.47)  or  (3.64),  we  have  to 
evaluate  integrals  of  the  form 


The  numerical  evaluation  of  the  integral  3 in  (4.1)  is  accomplished  by  choosing  a quad- 
rature formula  of  the  form 

m 

3=1  = Lf(xi)wi  (4.2) 

hi 

where  xi  are  the  integration  points  in  the  interval  X\  < x <x2 ; the  wi  are  numbers,  called 
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Figure  4.2:  Flow  Charts  of  Components  of  CODE 
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SETINT  = Calculate  quadrature  rule  data 

PREP  = Preprocessor:  Calls  routines  to  read  and  generate  data 

PROC  = Processor:  Forms  and  solves  finite  element  equations 

RCON  = Reads  control  parameters 

RNODE  = Reads  and  generates  nodal  point  coordinates 

RELEM  = Reads  and  generates  element  data 

RBC  = Reads  boundary  condition  data 

RMAT  = Reads  material  property  data 

RCONS  = Reads  input  constants  and  velocities 

RPARA  = Reads  parameters  for  Newmark  integration 

FORMGM  = Forms  global  coefficient  matrices 
MODGMB  = Adds  concentrated  mass  and  damp  at  top 
FORMGV  = Forms  global  forcing  vectors 

FMODGV  = Modifies  global  vectors  for  initial  prescribed  value  at  top  and  bottom 


APLEFO  = Adds  point  external  forcing 
APLBCV  = Eliminates  unnecessary  forcing  from  global  vector 
APLBCM  = Modifies  global  matrices  for  fixed  boundary  values 
NEWMRK  = Performs  Newmark  method  for  solving  time  series  response. 


See  chapter  6 for  details  about  algorithm,  iterations,  and  predictor-corrector 
scheme  using  extensibility  condition  for  vertical  degrees. 


ELEMM  — Calculates  elemental  coefficient  matrices 
ASSMBM  = Adds  element  matrices  to  global  matrices 

ELEMV  = Calculates  elemental  forcing  vectors 
ASSMBV  = Adds  element  vectors  to  global  vectors 

GETMAT  = Calculate  material  properties 
GETCON  = Gets  input  constants  and  parameters 
SHAPE  = Calculates  shape  function  values 

DRCHLT  = Imposes  boundary  conditions 

SOLVE  = Solves  the  system  of  linear  matrix  equation  for  unknown  displacement  field 


POST  = Prints  the  required  results 
TRI  = Performs  Gaussian  elimination 

RHSUB  = Performs  elimination  and  back  substitution  on  right  hand  side  of  equation 


boundary 


see  chapter  6. 


Figure  4.3:  A Brief  Description  of  Subroutines 
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CONTROL  PARAMETERS 

COMMON/CCON/NNODE,NELEM,NMAT,NCONS,NPOINT,NBW,NF 

NNODE  = number  of  node 

NELEM  = number  of  element 

NMAT  = number  of  materials 

NCONS  = number  of  characteristic  constants 

NPOINT  = number  of  concentrated  forces 

NODAL  POINT  DATA 

COMMON/CNODE/Z(I),X(I) 

Z = coordinates  of  nodal  points 
X = displacements  and  rotations  in  horizontal  direction. 

ELEMENT  DATA 

COMMON/CELEM/KIND(I),MAT(I),NCON(I),NODES(J,I),NINT(I) 

KIND  = degree  of  polynomial  in  element  shape  function 

MAT  = number  of  material  type  of  element 

NCON  = number  of  constant  type  of  element 

NODES  = node  numbers  of  nodes  in  element 

NINT  = number  of  integration  points  in  element  calculations 

MATERIAL  TYPE 

COMMON/CMATL/PROP(I,J),CONS(I,J) 

PROP  = material  properties 
CONS  = input  constants  and  velocity 

BOUNDARY  CONDITION  DATA 

COMMON/KBC(J,I),VBC(J,I),KPOINT(I),POINT(J,I),NB(I,J) 

KBC  = boundary  condition  flag  at  each  end 
VBC  = values  of  boundary  condition  data 
KPOINT  = node  number  of  point  loads 
POINT  = value  of  point  loads, 

NB  = prescribed  value  flag;  internally  decided 


Figure  4.4:  Data  Arrays  in  CODE 
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weights;  and  N/  is  the  order  of  formula.  Again,  we  can  replace  these  calculations  as  be- 
ing performed  on  a master  element  O defined  by  the  interval  -1  < £ < 1 , so  that  limits 
of  integration  are  -1  and  1.  Specifically,  if  we  choose  as  the  transformation  from  E,  to  x 
the  relation 


then  we  have  * = when£  = -l,  x = x2when£  = 1,  and dx  = (dx/dQdZ,  = ±(x2 -xx)dEs. 
The  form  of  the  quadrature  formula  (4.2)  now  becomes 


The  values  of  £/,  and  w/  for  / = 1 N/  completely  determine  the  numerical  integra- 
tion rule.  In  this  problem,  Gaussian  quadrature  rules  are  used.  The  Gauss  rule  of  order  N/ 
will  integrate  exactly  polynomials  of  degree  2 N,  - l.The  integration  such  as  (4.1)  are 
performed  in  CODE  using  Gaussian  quadrature  of  order  1 ,2,3,  or  4.  The  data  are  stored 
in  common  block  CINT  by  subroutine  SETINT. 

SHAPE  which  is  called  by  ELEMM  and  ELEMV  calculates  the  values  of  the  shape 
functions,  their  first  and  second  derivatives  with  respect  to  the  master  element  coordinate 
at  a specified  value  of  £,.  This  routine  is  supplied  with  the  degree  of  freedom  in  plane  N 
and  the  value  of  E,  at  which  the  values  are  to  be  calculated.  The  calculated  values  are  re- 
turned in  arrays  PSI,  DPSI,  and  DDPSI. 

The  evaluation  of  the  integration  of  the  form  (4.4)  is  performed  in  ELEMM  and 
ELEMV.  The  values  of  the  material  properties  and  the  input  constants  which  contain  in 
the  integrands  are  supplied  by  GETMAT  and  GETCON.  Since  the  integrand  contains  the 
first  or  second  derivative  with  respect  to  the  global  coordinate  x,  we  use  the  chain  rule  to 
calculate  these  quantities  at  the  integration  points.  Thus, 


x = + l-(x2  + £) 


(4.3) 


(4.4) 


dNj  _ dNj_^=  dN  2 
dx  dt,  dx  dE,  x2~xi 


(4.5) 
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dr2  dt,2  ' dx  ' ^2  Vx2-*1 


(4.6) 


For  instance,  the  explicit  formulas  for  evaluating  the  integral  multiplying  the  se- 
cond derivative  of  shape  function  into  its  integrand  can  be  evaluated  as 


GETMAT  and  GETCON  routine  calculates  the  value  of  the  material  properties  and 
input  parameters  at  a given  point  respectively,  and  are  called  ELEMM  and  ELEMV. 


The  preprocessing  unit  of  CODE  exists  to  supply  the  data  required  for  the  calcula- 
tion of  the  coefficient  matrices  and  load  vectors  for  the  problem  to  be  solved.  The  defini- 
tion of  the  various  groups  of  data  is  made  in  subroutines  called  by  PREP. 

Routine  RCON  reads  the  value  of  the  control  parameters  NNODE,  NELEM, 
NMAT,  and  NPOINT.  (See  Fig.  4.4  for  definitions  of  these  parameters.)  The  allowable 
range  of  control  parameters  is  limited  by  the  storage  capacity  of  array  of  the  computer. 

Routine  RNODE  reads  data  that  allow  the  definition  of  the  coordinates  of  all  nodal 
points.  The  nodes  are  numbered  consecutively  from  1 through  NNODE.  The  coordinates 
are  stored  in  array  Z and  transmitted  and  nodes  are  generated  with  equal  spacing  in  vari- 
ous sections  of  the  grid. 

Routine  RELEM  reads  the  data  that  allow  the  definition  of  all  the  elements  in  the 
finite  element  model.  The  elements  are  numbered  consecutively  from  1 through  NELEM. 

The  data  that  define  element  number  I are  introduced  in  the  following  way  of  arranging: 

KIND(I)  = kind  of  element  (degree  of  polynomial  in  the  shape  function) 


9?  = 


(4.7) 


4.3  Preprocessing  Routines 


3 for  cubic  elements 
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MAT(I),  NCON(I)  NINT(I)  = See  Fig.  4.4. 

NODES(J,I)  J = 1,2,...,(KIND(I)+1)  = nodal  point  numbers  of  nodes  in  element 

These  data  are  stored  in  the  array  in  common  block  CELEM  and  transmitted  as  ne- 
eded through  this  block.  RELEM  should  read  data  records  that  allow  the  definition  of  the 
data  described  above.  It  is,  of  course,  not  desirable  to  require  a separate  data  record  to  be 
read  for  each  element  being  defined  same.  The  following  format  of  data  records  will  in- 
troduce the  complete  specification  of  element  data. 

RECORD  1 : NREC  = number  of  element  data  records 

RECORD  2:  KINDI,  MATI,  NCONI,  NINTI,  NODEI,  NUMBER 

KINDI  = value  of  KIND  for  elements  to  be  generated 

MATI  = value  of  MAT  for  elements  to  be  generated 

NCONI  = value  of  NCON  for  elements  to  be  generated 

NINTI  = order  of  the  integration  rule 

NODEI  = nodal  point  number  for  first  node  of  first  element 
to  be  generated 

NUMBER  = number  of  elements  to  be  generated  in  sequence; 

element  number  and  nodal  point  numbers  are  incremented 
as  elements  are  generated. 

The  coding  of  RELEM  must  cause  record  1 to  be  read,  then  a loop  entered  until 
record  2 is  read  and  the  appropriate  number  of  elements  generated. 

Routine  RBC  reads  data  that  defines  the  boundary  condition  at  each  end  of  the  re- 
gion. For  each  boundary  condition,  the  data  consists  of  a flag  KBC  and  values  VBC.  The 
meaning  of  these  variables  is  shown  below.  Note  that  I = 1 for  the  bottom  end  and  I = 2 
for  the  top  end. 

KBC(J,I)  = 1 for  boundary  with  the  prescribed  value,  otherwise  0 

VBC(J,I)  = The  prescribed  values  (J  = odd;  displacements,  even;  rotations) 
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Routine  RBC  also  reads  the  data  point  loads  (if  any).  For  each  of  the  NPOINT  loads,  the 
nodal  point  number  at  which  the  loads  acts,  KPOINT,  and  the  value  of  the  load,  POINT, 
are  to  be  read.  Further,  the  attached  mass,  damp,  and  spring  at  the  top  end  are  read  in 
both  lateral  and  rotational  degrees  if  it  is  required.  The  capital  letter  R within  the  vari- 
ables indicates  the  rotational  degree. 

The  descriptions  of  routine  RMAT,  RCONS,  and  RPARA  are  given  in  Fig.  4.3  and 
the  meaning  of  the  variables  is  written  below. 

PROP(J,I)  = material  properties;  I = element  number 

J = 1 ; mass  per  unit  length  (riser  mass  + added  mass  + mud  mass) 

2;  mass  of  mud  per  unit  length 
3;  density  of  sea  water 
4;  section  area  of  riser 
5;  outside  diameter  of  riser 
6;  the  value  of  top  tension 
7;  the  value  of  bottom  tension 
8;  coefficient  of  elasticity  E 
9;  moment  of  inertia  I 

CONS(J,I)  = input  constants  and  velocities;  I = element  number 

J = 1 ; the  velocity  of  internal  fluid 
2;  current  velocity  (if  necessary) 

3;  Cd 
4;  Ca 

need  more  constants  for  the  analysis  of  vortex-induced  vibration 
or  vibration  due  to  wave  or  current  loading  (See  section  6.3  or  6.4) 


DT  = time  interval  for  integration 
TMAX  = maximum  time  to  be  considered 
TF  = ending  time  for  point  force  (if  any) 

ALPHA  = parameter  value  for  Newmark  integration 
DELTA  = parameter  value  for  Newmark  integration 
ER  = error  ratio  for  convergence  = (new  value-  old  value)/old  value 
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4.4  Finite  Element  Calculation  Routines 

This  unit  consists  of  the  calling  routine  PROC  and  the  routines  that  calculates,  as- 
semble, modify  for  boundary  conditions,  and  solve  the  system  of  matrix  O.D.Es 
(ordinary  differential  equations).  PROC  contains  the  calls  to  the  various  components  of 
the  analysis  procedure,  which,  in  CODE,  consists  of  successive  calls  to  routines 
FORMGM,  MODGMB,  FORMGV,  FMODGV,  APLEFO,  APLBCV,  APLBCM,  and, 
NEWMRK. 

Routine  FORMGM  or  FORMGV  directs  the  calculation  of  element  coefficient  ma- 
trices or  load  vectors  evaluated  by  routine  ELEMM  or  ELEMV  and  the  accumulation 
into  the  global  matrix  or  vector  by  ASSMBM  or  ASSMBV.  These  routines  begin  by 
zeroing  the  global  matrices  and  load  vector  respectively.  Then,  element  by  element,  the 
routines  that  calculate  and  assemble  the  element  matrices  (EM,  EC,  EK),  and  vectors 
(EFC),  are  called.  The  assembly  routines  (ASSMBM,  ASSMBV)  require  the  element  and 
global  matrices  or  vectors,  and  the  nodal  point  numbers  (equation  numbers)  for  the  nodes 
in  the  element.  Since  the  node  numbers  for  each  element  are  stored  consecutively  in  the 
array  NODES,  it  suffices  to  transmit,  through  the  calling  sequence,  only  the  location  of 
the  first  node  number  of  the  element.  Then,  the  contributions  to  the  global  matrices  or 
vector  from  a single  element  are  added  to  the  accumulated  contributions  from  other  ele- 
ments. Specially,  in  the  accumulation  of  global  matrices,  full  bandwidth  technique  is  im- 
plemented to  reduce  the  size  of  global  matrices.  In  our  problem,  half-bandwidth 
technique  is  no  longer  adequate  because  some  of  our  element  matrices  is  not  symmetric. 

Routine  MODGMB  adds  the  concentrated  mass,  damp,  and  spring  at  the  top  end, 
by  the  accumulation  of  their  contributions,  to  the  corresponding  location  of  global  ma- 
trices and  APLEFO  adds  the  point  loads,  by  the  treatment  of  array  KPOINT  and 
NODES,  to  the  selected  position  of  global  vector.  In  some  type  of  boundary  condition 
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(no  movement  at  the  top  end  from  original  position),  the  operation  of  MODGMB  has  no 
effect  on  the  solutions  because  the  rows  and  columns  corresponding  to  top  boundary  are 
eliminated  in  subroutine  APLBCM.  Routine  APLBCM  and  APLBCV  eliminates  unnec- 
essary rows  and  columns  due  to  fixed  boundary  values  from  global  matrices  and  vector 
respectively. 

Routine  FMODGV  modifies  global  vectors  for  initial  prescribed  values  at  top  and 
bottom  boundary.  The  modification  for  boundary  conditions  is  made  by  calling  routine 
DRCHLT  representing  the  computer  implementation  of  the  essential  boundary  condi- 
tions, which  may  be  found  in  several  references  related  to  finite  element  book. 

Finally,  routine  NEWMRK  performs  Newmark  integration,  predictor-corrector 
scheme,  and  convergence  in  solution,  and  then  prints  results.  This  routine  will  be  dis- 
cussed in  chapter  6 in  details. 


CHAPTER  5 

SOLUTION  OF  THE  FREE  VIBRATION  OF  THE  SYSTEM 


For  the  investigation  of  the  characteristics  of  a marine  riser  system  with  the  inclu- 
sion of  internal  flow,  it  is  generally  required  to  estimate  the  system  frequencies  and  mode 
shapes.  However,  it  is  not  a simple  task  for  the  case  of  the  system  including  highly  non- 
linear terms.  Moreover,  even  the  attempt  to  resolve  the  problem  has  not  been  made  for 
the  case  of  the  system  with  deformation-dependent  nonlinear  terms  like  the  present  case. 
Thus,  the  system  frequencies  and  mode  shapes  are  obtained  by  solving  the  approximate 
form  of  governing  equations  derived  in  section  2.5.3  or  3.5.  This  serves  our  purpose  to 
investigate  the  effect  of  internal  flow  on  the  system. 

This  chapter  consists  of  three  sections;  first,  semi-analytical  solutions  are  obtained 
using  series  expansion  for  the  mathematical  model  and  second,  numerical  solutions  are 
calculated  adapting  complex  eigenvalue  method  for  the  matrix  equilibrium  equation.  Fi- 
nally, the  algorithm  for  the  computer  program  is  introduced  and  verification  is  given  by 
comparing  the  solutions  from  the  first  method  to  that  from  the  second  method.  Further, 
the  comparisons  of  semi-analytical  solutions  with  results  from  other  papers  is  presented. 

5.1  Semi-Analvtical  Solution 

The  approximate  governing  equation  for  lateral  vibration  of  a marine  riser  system 
with  the  internal  flow  and  linear  tension  has  been  derived  in  section  2.5.3.  The  free 
vibration  equation  of  the  system  can  be  derived  simply  by  setting  the  forcing  term  in  the 
right  hand  side  of  Eq.  (2.68)  to  zero  and  dropping  the  subscript  1 of  variable  x indicating 
the  lateral  deformation.  The  resulting  equation  represents  the  free  vibration  of  the  system 
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in  the  lateral  direction  and  is 

m,x  + 2mf\ix'  + m/V \x"  + (EIx")"  - wx'  - Tex"  =0  ( 5.1  ) 

where  the  unknown  variable  x represents  the  displacement  in  lateral  direction  and  prime 
denotes  the  derivative  with  respect  to  a vertical  coordinate  z. 

Providing  that  pipe  has  an  uniform  section  of  constant  El  and  selecting  Te  as  the 
mean  value  of  the  tension  force  acting  along  the  length  of  riser,  Eq.  (5.1)  becomes 

m,x  + 2m /V Ix/  + mfV]x"  + EIx""  - wx'  - Tex"  = 0 ( 5.2  ) 

where  f e represents  mean  tension  and  is  given  as  TTB  + wL/2  or  TTR  - wL/2. 

5.1.1  Quasi  Mode  Shape 

Equation  (5.2)  differs  from  the  usual  vibration  equations  in  that  it  contains  a mixed 
derivative  term  with  respect  to  time  and  space  which  is  the  coriolis  force,  the  force  re- 
quired to  rotate  fluid  elements  with  local  pipe  rotation.  Its  mixed  derivative  causes  an 
asymmetric  distortion  of  classical  mode  shapes  and  could  lead  to  a flutterlike  instability. 
Thus,  Eq.  (5.2)  does  not  possess  the  classical  normal  modes.  Its  solution  can  not  be  sepa- 
rated simply  into  time  and  space  components.  For  example,  if  a trial  solution  of  the  form 

x(z,i)  = x(z)  sin©/  (5.3) 

is  substituted  into  (5.2),  it  can  be  seen  that  the  coriolis  force  term  varies  as  cosofr  , while 
the  remainder  of  the  terms  have  a sinco?  term.  This  suggests  that  solution  should  be  writ- 
ten as 

x(z, /)  = aix(z)  sin  q/  + a2x  (z)  cos  ofr  ( 5.4  ) 

The  boundary  conditions  for  a pinned-pinned  span  are  given  by 

x(0, t ) = x(L,/)  = d2x(0 ,t)ldz2  = 52x(L ,t)/dz2  = 0 ( 5.5  ) 

and  those  boundary  conditions  are  satisfied  by  the  set  of  sinusoidal  mode  shapes, 


x (z)  = sin  rmz/ L, 


«=  1,2,3, ... 


(5.6) 
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These  mode  shapes  make  the  first,  third,  fourth,  and  sixth  terms  in  Eq.  (5.2)  unaltered, 
but  the  mixed  derivative  in  the  coriolis  force  term  generates  spatially  asymmetric  terms 
for  a symmetric  mode  shape  (n  = 1,3, 5,  • • •)  and  spatially  symmetric  terms  for  an  asym- 
metric mode  shape  ( n = 2,4, 6,  • • •).  These  considerations  imply  that  the  solution  of  (5.2) 
with  pinned  end  conditions  should  be  the  sum  of  symmetric  and  asymmetric  spatial 
modes  with  sine  and  cosine  components  (Housner,  1952;  Blevins,  1990). 

Xj(z,0  = X a„  sin  («7tz/L)sincoj/ + X a„  sin  (wtz/L)coscoji 

"=1,3,5,- ••  «=2,4,6,- 

j = 1,2,3,--  (5.7) 

where  coj  is  the  natural  frequency  of  the  j-th  vibration  mode  and  the  coefficient  equation 
will  be  derived  in  next  section. 

Eq.  (5.7)  shows  that  the  classical  mode  shape  does  not  exist  in  our  system  as  stated 
early.  Moe  (1988)  has  shown  that  there  is  a phase  shift  due  to  the  internally  flowing  flu- 
id. Chen  has  made  a heuristic  assumption  that  the  time  components  have  been  dropped 
out  of  equation  (5.7)  at  the  end  of  his  derivation  and  showed  through  his  numerical  ex- 
ample that  the  error  due  to  the  assumption  was  negligible  as  compared  to  other  author's 
results  from  different  approaches.  According  to  his  assumption,  the  1 st  quasi  mode  shape 
can  be  formed  by  dropping  the  time  components  from  (5.7)  and  given  by 

*i  (z)  - X an  sin  (wiz/L)  + X a„  sin  (wtz/L)  ( 5.8 ) 

n=l,3,5,-  »=2,4,6,-  • • 

where  the  coefficients  was  determined  to  the  1st  natural  frequency  co  i. 

In  this  paper,  the  quasi  mode  shape,  which  is  not  from  Chen's  assumption  but  by 
including  time  dependent  component  in  Eq.  (5.7),  is  implemented  for  the  determination 
of  system  frequency. 
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5.1.2  Modified  Eigenvalue  Problem 

In  order  to  obtain  the  coefficient  equations,  let's  begin  with  the  substitution  of  trial 
solution  (5.7)  into  Eq.  (5.2).  By  this  substitution,  the  coriolis  force  and  the  fifth  term  in 
(5.2)  produce  terms  containing  cos  (nnz/L) . These  cosine  terms  can  be  expanded  in  a 
Fourier  half-range  series  of  sine  functions 

cos(wuz/L)  = X bnp  sin  (pnz/L),  n=  1,2,3, •••  (5.9) 

p=  1,2,3,- •• 

where 


/ 0 

\ 4p/{n(p2  - n2)} 


n + p - even 
n + p = odd 


(5.10) 


This  Fourier  series  converges  relatively  slowly  over  the  span  and  not  at  all  at  the  ends, 
but  it  does  allow  the  spatial  dependence  to  be  factored  out  of  the  solution.  With  the  sub- 
stitution of  the  series  expansion  of  cosine  function,  the  terms  in  (5.2)  can  be  grouped  ac- 
cording to  whether  they  contain  sin  cot  or  cos  &t  to  give  following  equation: 


Z a„  -m/Vj(wt/L)2  + Te(wr/L)2  +EI(wr/L)4}sin(wuz/L)sin(©,7) 

«=1,3,5,-  • • v ' 

+ Z (-Sm/Vidj/L)  j Z (appn)/(n2 -p2)lsin(nnz/L)sin(Qjt) 


+ Z 

n=  1,3,5,- 


-w(«7t/L) 


I 

V=2,4,6,- 


Ap/(n(p2  - «2))sin  (pnz/L)  )sin(07/) 


+ Z a n {-m tQ j -ffifW i (nn/L)2  +Te(nn/L)2  +EI(nn/L)4} sin (nnz/L)cos((Ojt) 

n= 2,4,6,-  • • 


+ Z (SnifViQj/L)  j Z (appn)/(n2 -p2)lsin(nnz/L)cos(G)jt) 


n= 2,4,6,- ■■ 


(_P=1,3,5,--- 


+ Z a„{-w(«7t/L)  X 4p/(n(p2  - «2))sin (pnz/L)  }cos(co,/)  = 0 (5.11) 

»=2,4,6,---  ( Vp=l,3,5,---  J \ 
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For  a neutrally  buoyant  case  (w  = 0),  Eq.  (5. 1 1)  becomes 
^ <Ti  -/«/Vi(wr/L)2  +Te(wt/L)2  + EI(/m/L)4  }sin  (m^/L)sin (co//) 


+ X (— 8/My V !© j/L)  j X (appri)/(n2  -p2)  [-sin  (mcz/L)sin  (©.-/) 

”=1,3,5,- ••  [P=2,4,6,---  J 

X an  {-/M/CO2  -/m/V2(m7t/L)2  +Te(«7t/L)2  + EI(«7t/L)4}sin(«7tz/T)cos(©/) 


”=2,4,6,- 


+ X (8/m/Vjq/L))  X (appn)/(n2  -p2)  isin  (//7u/L)cos  (co ,/)  = 0 (5.12) 

”=2,4,6,-  ••  [/>=1,3,5,---  J 

The  coefficients  of  each  group  are  set  to  zero  to  give  the  following  equations 

a„[EI(/ni/L)  -/M/V2(«7i:/L)2 + Te(//7t/L) -/m,co2  J = (8/M/ViOy/L)  X appn/(n2  ~p2) 

p=2,4,6,--- 

az  = 1,3,5,  (5.13) 

^«[eI(/2ti/L)  - nif  V 2 (mi/L)2  + Te(/m/L)  - mt<x>j  J = -(8/m/V [Co^/L)  X appn!(n2  - p2) 


p=  1,3,5,- 


//  = 2,4,6,--- 


(5.14) 


For  a negatively  or  positively  buoyant  case  (w  * 0),  the  exact  coefficient  equations  can 
not  be  obtained  because  the  third  and  sixth  terms  in  (5.1 1)  have  different  running  num- 
ber in  summation.  It  causes  a slight  shift  in  phase  to  remaining  terms.  However,  its  mag- 
nitude is  small  enough  to  ignore  phase  lag  if  a large  number  of  summation  is  taken.  Thus, 
the  running  number  n and  p can  be  switched  in  the  summation  of  the  third  and  sixth 
term  of  (5.1 1)  with  each  other  and  then,  Eq.  (5.1 1)  becomes 


X an  {-/M/(o2  -/M/V2(/27t/L)2  + Te(«7r/L)2  + EI(//7t/L)4]sin (//7Cz/L)sin (®jt) 


+ X (-8/m/Vi© /L)  5 X (appn)/(n2  -p2)  fsin(//7tz/L)sin(o,A 

”=1,3,5,-  ••  (p=2,4,6,---  J 

f \ 

+ X -4w/L  X appn/(n2  -p2)  sin  (mtz/L)sin(©^) 

n=  1,3,5,*  ••  \p= 2,4,6,  ••  J 

+ X a„  {-/m,co2  - mz/ V 2 (mtc/L) 2 + Te(//7t/L)2  +EI(//7t/L)4}sin(/?7tz/L)cos(ro77) 

”=2,4,6,-  ■ • 
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+ Z (8/w/Vi©/L)  j Z (appn)/(n2  -p2)  is  in  (wtz/L)cos  (©/) 


n=  2,4,6, 


lp=  1,3,5, 


r \ 

+ Z -4w/L  Z appn/(n2  -p2)  sin(wtz/L)cos(©/) 

" " ' V=i,3,5,-  • • y 


"=2,4,6,- 


(5.15) 


Collecting  Eq.  (5.15)  into  groups  according  to  sin©/ or  cos©/ and  setting  the  coeffi- 
cient to  zero,  the  following  set  of  coefficient  equations  can  be  obtained: 

a„[EI(mt/L)4  -mfV\(rmfL)2  +Te(wr/L)  -mtG>j  ] 

= (8  mfV  j (Oj  + 4w ) Z apPnl(n 2 - p 2)  «=1,3,5,---  (5.16) 

a„[EI(wt/L)4  -/w/V2(wt/L)2  + Te(wt/L)  - mtQj  ] 

= (-8/w/Viffly  + 4w)  Z appn/(n2-p 2)  w = 2,4,6,---  (5.17) 

p=l,3,5,-  •• 

These  equations  can  be  put  in  matrix  form: 

{[^np]  — ©y  [C„p]  - [/]]  (a } = 0 (5.18) 

where, 


K, 


np 


EI(rat/L)  4 + T e (nn/L) 2 -m/V2  (wt/L)  2 , 
(-4 w/L)  {/?«/(« 2 -p 2 ) } , 
(-4w/L) {pn/(ri 2 -p2)}. 


n-p 

n = odd,  p = even 
n = even,  p = odd 
n^p,  «+/?  = eve  n 


(5.19) 


0,  n -p 

Cnp  - (&mfVifL){pn/(n2  -p2)},  « = odd,  /?  = even 

(-8/w/Vi/L){p«/(«2  -p2)},  « = even,  p = odd 
v 0,  n*p,  n+p  = even 

(a)  = {ai,a2,a3, }T 


[I  ] = the  identity  matrix  with  value  of  one  on  the  diagonal  and  all  other  equal  to  zero. 
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Nontrivial  solution  of  Eq.  (5.18)  are  sought  by  setting  the  determinant  of  the  coef- 
ficient matrix  to  zero: 

| [£]  - © [C]  - co2m,[/]  | = 0 (5.21) 

Because  the  system  has  an  infinite  number  of  natural  modes,  we  have  to  consider  only 
the  first  few  modes  practically.  Blevins  (1990)  included  only  the  first  two  modes  in  his 
appropriate  analysis,  then  #3,  a4;  a$  - ■ ■ are  set  equal  to  zero.  In  this  paper,  the  more  num- 
ber of  natural  modes  will  be  included  in  practical  analysis. 

5.2  MDOF  System  with  Coriolis  Matrix 

Since  the  semi-analytical  method  included  the  assumption  in  derivation  and  the 
simplicity  in  a practical  sense  as  shown  in  the  previous  section,  the  error  resulting  from 
the  assumption  and  the  simplicity  should  be  checked.  The  conventional  way  for  deter- 
mining the  system  frequency  is  to  construct  eigenvalue  problem  using  the  matrix  equilib- 
rium equation  of  free  vibration  system  which  is  already  constructed  in  section  3.5  as 
given  the  following  equation: 

[M] {*(/)}  + [Cj ] (x(0 } + [Tij - Cej  + kj]{x(t)}=  0 (5.22) 

Rewriting  (5.22)  in  a usual  simple  form,  we  get 

MU+CU+KU=0  (5.23) 

In  above  equation,  C is  not  a usual  damping  matrix  but  a coriolis  matrix  which  is  skew- 
symmetric.  As  discussed  earlier,  C U acts  like  a dynamic  coupling  force  not  a damping 
force.  Thus,  it  can  not  be  allowed  to  construct  eigenvalue  problem  simply  using  Duncan's 
method  (1956)  because  the  imaginary  part  of  eigenvalue,  which  indicates  a damped  fre- 
quency, is  always  zero.  The  complex  value  with  zero  imaginary  should  be  implemented 
and  then,  the  imaginary  part  has  to  be  checked  whether  it  is  equal  to  zero.  This  is  verified 
through  the  numerical  example  given  in  5.3.2. 
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For  the  construction  of  eigenvalue  problem  using  Duncan's  method,  we  combine 
the  identity 


MU-MU  = 0 
with  Eq.  (5.23)  to  get 


0 M 
M C 


U 

u 


> + 


-M  0 
0 K 


U 

u 


0 

0 


Equation  (5.25)  can  now  be  rewritten  as 

79t7(,  + X7l  = 0 

where 

7ft  = 


0 M 
M C 


7t 


-M  0 
0 K 

U 

u 


I 

u 

r and 

7C  = 

> 

u 

Now  letting 
71  = \ept 

where  p is  complex.  Eq.  (5.26  ) can  be  transformed  into 
( j>7tt  + X)\  = 0 


(5.24) 


( 5.25  ) 


(5.26) 


(5.27) 


(5.28) 


(5.29) 


This  last  equation  is  of  a standard  form  for  which  eigenvalue  computer  programs 
are  available.  The  disadvantage  we  have  for  this  simplification  is  that  the  sizes  of  all  ma- 
trices are  doubled. 
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5.3  Numerical  Evaluation  of  System  Frequency 
5.3.1  Algorithm  for  Computer  Program 

In  the  previous  section,  eigenvalue  problem  with  the  form  of  (5.18)  and  (5.29)  has 
been  derived.  Eq.  (5.29)  is  of  a standard  form  for  which  the  algorithms  to  obtain  the  re- 
quired eigenvalues  and  eigenvectors  are  available,  on  the  other  hand,  (5.18)  is  not  in  the 
standard  form.  Thus,  it  is  required  to  implement  a slightly  different  algorithm  based  on 
one  for  the  standard  form  of  eigenvalue  problem.  As  recognized  in  the  algorithm  summa- 
rized in  Fig.  5.1,  the  standard  form  of  eigenvalue  problem  has  to  be  resolved  for  every 
iteration  step. 


solve  { [A](i)+^  [B] } {U}  = 0 
obtain  co(iit^1/2 


NO 


Obtain  eigenvalue  and 
corresponding  eigenvector 


Let 

oi  * 5 = o)ii) 


Figure  5.1:  Algorithm  for  Solving  Non-Standard  Eigenvalue  Problem 
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In  this  paper,  the  inverse  iteration  algorithm  with  shifting  value  is  chosen  to  solve 
the  given  eigensystem.  The  technique  of  inverse  iteration  is  very  effectively  used  to  cal- 
culate an  eigenvector,  and  at  the  same  time  the  corresponding  eigenvalue  can  also  be 
evaluated.  Inverse  iteration,  like  a usual  vector  iteration  method,  can  be  employed  even  if 
some  diagonal  elements  of  mass  matrix  has  zero  value  as  shown  in  Eq.  (5.27).  Further, 
the  shift  of  eigenvalue  improves  the  convergence  rate  in  the  vector  iteration  and  acceler- 
ates the  calculation  of  the  eigenvectors.  Those  techniques  are  utilized  together  in  this  pa- 
per and  the  computer  program,  of  which  algorithms  are  referenced  from  the  text  of  Bathe 
and  Wilson  (1976),  is  developed  in  this  study.  In  addition  to  the  subroutines  in  computer 
program,  subroutine  FORMA  which  constructs  the  coefficient  matrices  in  (5.18)  and  first 
step  of  Fig.  5.1,  and  FORMB  which  forms  the  modified  stiffness  matrix  [ A ](,)  for  every 
iteration  step  as  shown  in  Fig.  5.1  are  built  in  computer  program.  While  FORMA  and 
FORMB  are  necessary  routines  for  the  calculation  of  the  eigensystem  by  first  method 
(semi-analytical  method),  FORMA  is  only  routine  for  preparing  coefficient  matrices  Ttt 
and  7:  in  Eq.  (5.29),  whose  elements  are  filled  with  the  resulting  coefficients  from  the 
matrix  equilibrium  equation  (Eq.  (5.22)  and  (5.23))  constructed  by  finite  element  approx- 
imation (see  section  3.5). 

5.3.2  Comparison  of  Results 

For  the  verification  of  the  solutions  from  semi-analytical  method,  the  comparisons 
are  made  in  two  different  ways.  The  first  is  to  calculate  the  system  frequencies  of  the  sys- 
tem with  the  exclusion  of  internal  flow  using  semi-analytical  method  and  then  compare 
those  with  the  results  from  previous  investigators.  This  comparison  provides  the  validity 
of  the  representation  of  solution  into  the  sum  of  two  components  as  shown  in  Eq.  (5.7) 
but  it  does  not  supply  enough  evidence  for  the  system  including  the  internal  flow.  Thus, 
another  one  in  section  5.2  is  implemented  to  compare  with  the  results  from  the  semi- 
analytical  method.  This  comparison  gives  enough  evidence  even  for  the  system  including 
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internal  flow.  The  comparison  with  Chen’s  results  ( 1 992)  may  be  also  given.  However,  in 
the  expansion  of  cosine  term  into  Fourier  series,  he  followed  the  wrong  use  of  Fourier 
coefficients  which  was  also  found  in  Housner's  paper  (1952).  The  wrong  expansion  of 
the  coriolis  force  term  may  lead  to  the  error  in  the  calculation  of  the  system  frequency 
with  internal  flow. 

The  example  data,  which  is  for  a typical  drilling  riser  system  for  use  in  the  northern 
North  sea  and  repeated  from  Chen  (1992),  is  given  in  Fig.  5.2.  Using  those  data,  the  first 
comparison  is  given  in  Table  5.1  which  presents  the  comparison  with  the  natural 
frequencies  of  the  1st  to  10th  modes  obtained  from  semi-analytical  method,  Dareing 
(1976),  Spanos  (1980),  Kim  (1988),  and  Chen  (1992).  The  discrepancy  between  those 
methods  including  semi-analytical  method  is  negligible,  that  is,  the  representation  of 
solution  such  as  Eq.  (5.7)  is  acceptable. 


Outside  diameter 
Constant  of  elasticity 
Sectional  moment  of  inertia 
Riser  length 
Riser  mass 


D - 2.0  ft,  with  5/8"  wall  thickness 
E = 30  x 106  psi 
I = 0.1518  ft4 
L = 500  ft 
m = 20.86  slug/ft 

(includes  mass  of  drilling  mud  and  sea  water) 


Bottom  tension  TTB  = 

Effective  weight  per  unit  length  w = 
Mean  tension  Te  = 

Density  of  drilling  mud  pm  = 

Density  of  sea  water  pw  = 


286,000  lb 
264.68  lb/ft 
326,830.5  lb 
85  lb/ft3 
64.8  lb/ft3 


Figure  5.2:  The  Design  Properties  and  Data  for  a Drilling  Riser  System 
for  Use  in  the  Northern  North  Sea 
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Table  5.1:  Comparison  of  System  Frequencies  of  the  1st  to  10  th  Modes  for  a Drilling 
Riser  System  with  No  Internal  Flow.  (System  frequency  = co;  (rad/sec),  Vi  = 0 m/sec) 


Mode 

i 

Semi-Analytical 

Method 

Chen 

Dareing 

Kim 

Spanos 

1 

0.81813 

0.81585 

0.81498 

0.81484 

0.831 

2 

1.80520 

1.80430 

1.80362 

1.80390 

1.831 

3 

3.08798 

3.08585 

3.08762 

3.08328 

3.123 

4 

4.73703 

4.73268 

4.73748 

4.73328 

4.778 

5 

6.78865 

6.78123 

6.78901 

6.78934 

6.832 

6 

9.26100 

9.24969 

9.301 

7 

12.1634 

12.14743 

12.193 

8 

15.5008 

15.47942 

15.506 

9 

19.2760 

19.24850 

19.232 

10 

23.4907 

23.45631 

23.352 

For  the  riser  system  with  internal  flow  of  6 m/sec,  the  comparisons  of  the  1st  to 
10th  modal  frequencies  calculated  using  both  semi-analytical  method  and  numerical 
method  described  in  section  5.2  are  given  in  Table  5.2.  The  system  frequencies  obtained 
using  second  method  are  complex  values  which  consist  of  real  and  imaginary  part.  The 
imaginary  part  represents  damping  frequency.  As  shown  in  Table  5.2,  the  imaginary 
parts  have  almost  zero  values,  that  is,  the  mixed  derivative  term  in  governing  equation  or 
the  coriolis  skew- symmetric  term  in  matrix  equilibrium  equation  does  not  have  the  same 
effect  as  that  of  viscous  damping  force  although  it  involves  the  time  derivatives.  In  other 
words,  the  mixed  derivative  term  or  the  coriolis  matrix  indicates  the  dynamic  coupling 
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Table  5.2:  Comparison  of  System  Frequencies  of  the  1st  to  10  th  Modes  for  a Drilling 
Riser  System  with  Internal  Flow.  (System  frequency  = ©i  (rad/sec),  Vt  = 6 m/sec) 


Mode 

i 

Semi-Analytical 

Method 

Numerical  method  described  in  section  5.2 

Real  value 

Imaginary  value 

1 

0.81701 

0.81417 

0.00231 

2 

1.80326 

1.80030 

0.00194 

3 

3.08550 

3.06014 

0.00102 

4 

4.73421 

4.73527 

0.00032 

5 

6.78562 

6.78776 

0.00003 

6 

9.25782 

9.26105 

0.000007 

7 

12.1601 

12.1916 

0.000001 

8 

15.4975 

15.4908 

0.000000 

9 

19.2727 

19.2735 

0.000000 

10 

23.4889 

23.4909 

0.000000 

force  as  discussed  earlier.  Further,  it  can  be  recognized  from  the  table  that  the  discrepan- 
cy in  the  system  frequencies  resulting  from  the  two  different  methods  is  negligible.  Thus, 
both  two  methods  developed  in  this  chapter  are  acceptable  for  the  estimation  of  the  sys- 
tem frequency  regardless  of  the  existence  of  internal  flowing  fluid. 

Of  two  methods,  semi-analytical  method  is  more  convenient  than  the  other  because 
there  is  no  need  to  construct  the  coefficient  matrices  through  the  finite  element  approx- 
imation. Thus,  semi-analytical  method  will  be  adapted  for  the  investigation  of  the  effect 
of  internal  flow  on  system  frequency. 


CHAPTER  6 

FORCED  VIBRATION  OF  THE  SYSTEM 


This  chapter  deals  with  the  forced  vibration  of  the  system  derived  in  chapter  3 and 
its  dynamic  behaviors.  The  system  considered  here  consists  of  two;  first,  the  system  sim- 
ply responding  to  the  general  environmental  loads  such  as  wave  or  current  and  second, 
the  system  responding  to  the  vortex  loads  generated  in  the  transverse  direction  to  inline 
current,  i.e.,  vortex-induced  vibration.  Newmark  method  is  used  for  time  step  integration 
and  predictor-corrector  scheme  is  implemented  as  an  iteration  algorithm  for  the  nonlinear 
terms.  Using  the  techniques  mentioned  above,  vortex-induced  vibration  due  to  current 
and  inline  vibrations  due  to  current  or  wave  are  obtained.  And  also,  the  application  of  top 
boundary  condition  to  the  system  matrices  or  vectors  is  discussed  and  its  effect  to  the  re- 
sponses of  a riser  is  inspected. 

6.1  Newmark  Method  for  Time  Step  Integration 
In  chapter  3,  we  derived  the  equations  of  motion  for  a riser,  which  can  be  written 
in  the  following  general  form 

MU  + CU  + KU  = R ( 6.1  ) 

where  M,  C,  and  K are  coefficient  matrices;  R is  external  load  vector;  and  U,  U and  U 
are  the  displacement,  velocity,  and  acceleration  vectors  of  the  finite  element  assemblage. 
However,  it  should  be  noted  that  although  Eq.  (6.1)  seems  to  be  linear  in  its  representa- 
tion, it  is  actually  nonlinear  because  the  equation  is  written  in  the  abridged  form  of  Eq. 
(3.51),  (3.61),  (3.62),  and  (3.63).  Thus,  the  coefficient  matrices  or  vectors  in  (6.1)  are 
also  determined  by  comparing  (6.1)  with  (3.51),  (3.61),  and  so  on. 
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Previous  studies  (Bathe  and  Wilson,  1976)  have  shown  that  Newmark  method  is 
the  most  convenient  and  economical  one  to  employ  a time  step  integration  because  larger 
time  steps  can  be  used  and  still  obtain  accurate  values  for  the  response.  Its  stability  analy- 
sis has  been  done  for  linear  systems.  They  do  not,  however,  carry  over  to  nonlinear  case 
up  to  now.  Unfortunately,  a relatively  large  time  step  may  lead  to  the  divergence  in  solu- 
tion. Thus,  it  is  generally  suggested  to  follow  the  Newmark’s  algorithm  developed  for 
linear  case  and  apply  simple  iteration  outside  the  loop  of  its  flow  or  another  scheme  for 
iteration  (See  next  section.)  to  remove  the  possibility  of  unbounded  solution. 

The  Newmark  integration  scheme  is  an  extension  of  the  linear  acceleration  method. 
The  following  assumptions  are  used: 


where  a and  bare  parameters  that  can  be  determined  to  obtain  integration  accuracy  and 
stability.  When  5 = \ and  a = ±,  relations  (6.2)  and  (6.3)  correspond  to  the  linear  accel- 
eration method.  Newmark  originally  proposed,  as  an  unconditionally  stable  scheme,  the 
constant  average  acceleration  method,  in  which  case  6 = \ and  a = \ (See  Fig.  6.1). 


(6.2) 


(6.3) 


Figure  6.1:  Newmark's  Constant  Average  Acceleration  Scheme 
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In  addition  to  (6.1)  and  (6.2),  for  solution  of  the  displacements,  velocities,  and  ac- 
celerations at  time  t + At,  the  equilibrium  equations  (6.1)  at  time  t + At  are  also  consi- 
dered: 

MU/+a/  + CIW  + KU/+4/  = R,+a/  ( 6.4) 

Solving  from  (6.3)  for  Uh-a/  in  terms  of  Uh-a/,  and  then  substituting  for  U,+A/  into  (6.2), 
we  obtain  equations  for  and  U/+a/>  each  in  terms  of  the  unknown  displacements 
U/+A/  only.  These  two  relations  for  U,+A,  and  U,+A/  are  substituted  into  (6.4)  to  solve  for 
U)+A/,  after  which,  using  (6.2)  and  (6.3),  U/+a/  and  U/+Af  can  also  be  calculated.  The  com- 
plete algorithm  using  Newmark  integration  scheme  is  given  in  Figure  6.2. 


A.  Initial  Calculations: 

1 . Form  K,  M,  and  C. 

2.  Initialize  U0,  U0,  and  U0. 

3.  Select  time  step  size  At,  parameters  a and  5,  and  calculate  integration  constants: 

5 > 0.50,  a > 0.25(0.5  + S)2, 

a0  = l/(aA t2),  a,  = 5/(aA t),  a2  = l/(aA/),  a3  = 0.5a  - 1, 
a4  = 5/a  - 1,  a5  = 0.5A/(5/a-2),  a6  = At(l-5),  a7  = 5 At 

4.  Form  effective  stiffness  matrix  K:K  = K + a0M  + aiC. 


B.  For  Each  Time  Step: 

1.  Calculate  effective  loads  at  time  t + At  : 

A 

R/+a/  = R/+a/  + M (a0U/  + 02I. Jt  + ct3  U/ ) + C(aiU,  + a 4 U/  + c/5 1 1/ ) 

2.  Solve  for  displacements  at  time  t + At  : 

A A 

KUffAi  = R,+a, 

3.  Calculate  accelerations  and  velocities  at  time  t + At: 

Ur+A/  — a o (U<+a/  — U/)  — ct2  U t — ^3  U| 

U t+Al  = V t + 06  v 1 + Clj  U/+AJ 


Figure  6.2:  Algorithm  for  Step-by-Step  Solution  Using  Newmark  Integration  Method. 
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6.2  Predictor-Corrector  Scheme 

The  utilization  of  the  kinematic  constraint  on  unit  tangent  vector  makes  it  feasible 
to  represent  the  vertical  degrees  of  freedom  in  terms  of  the  lateral  ones  as  stated  in  chap- 
ter 2.  This  does  lead  to  the  reduction  of  the  degrees  of  freedom  per  element  from  12  to  8 
and  remove  the  possibility  of  the  divergence  in  solutions  due  to  the  iteration  of  highly 
nonlinear  terms  in  the  vertical  degrees  of  freedom.  Moreover,  it  can  be  implemented  as 
an  another  algorithm  scheme  for  the  iteration  due  to  nonlinear  terms.  For  the  develop- 
ment of  an  algorithm,  so  called  predictor-corrector  scheme,  the  derivation  of  the  incre- 
mental form  of  the  extensibility  condition  is  necessary  because  the  equations  of 
equilibrium  for  a finite  element  system  in  motion  to  be  solved  is  of  the  incremental  form. 

Using  the  kinematic  constraint  on  unit  tangent  vector  of  riser,  i.e.  It  ■ t = r'  ■ r'  = 1 , 
we  have  the  following  relationship: 

x\ 2 + x22  + (x3  + 1/(1  + £»))2  = 1 ( 6.5  ) 

Applying  the  incremental  operator  A on  above  equation,  we  obtain 

Axg  = - (x[Ax[  + X2AX2HX3  + 1/(1  + e'))-1  + Ae,(l  + e<)~2  ( 6.6  ) 

Also,  we  have  Te  = F • t = F • r'  = Fix^  + F2*2  + F3(*3  + 1/(1  + e,))  ( 6.7  ) 

Applying  the  incremental  operator  on  (6.7),  we  get 

ATe  = AF  \x\  + AF2*2  + AF3(*3  + 1/(1  + 8,))  + FiAx^  + F2AXj 

+ F3A*3  - F3As,(1  +S,)"2  ( 6.8  ) 

From  Eq.  (2.16)  and  (2.17),  we  have 

Te  = EAs/  + PwAo^^w  — S — x 3 ) — pmAjg(Sm  — S — xj ) ( 6.9  ) 

Applying  A again  on  (6.9),  we  obtain 


ATe  = EA As,  - Ax3g(p„,A0  - pmA,) 


(6.10) 
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Eliminating  ATe  from  (6.9)  and  (6.10),  we  compute  As, 


AFix1+AF2X2+AF3(x3+l/(l+e,)j+FiAxj+F2Ax2+F3Ax3+Ax3g-(pM,Ao-pOTA/) 


As,  = __ 

EA  + F3  (1  +s,)-2 

Substituting  (6.1 1)  into  (6.6)  and  solving  for  Axj , we  obtain 

/ . / / / 
f3 


(6.11) 


Ax'  = (l 


+ 


x XjAxj+XjAxj 

/ " ~ZT 


+ 


1 


{ AFjx'  + AF2x'2 


EA(1+s*)2/  x3  + 1/(1  +s,)  EA(l+s;)2+F3 

+ AF3(x'  + 1/(1  + s,))  + FiAx^  + F2Ax2  + Ax3g(pwA0  - p,„A,)}]  ( 6.12  ) 

In  addition  to  (6.12),  we  also  have 


Ax3  = J Ax 3 ds\ 
x3  = (1  -x'2-x'2)1/2  - 1/(1  +8,) 


^3  = x3(0)  + X3  ds  j 


r 

Jo 


(6.13) 


As  stated  earlier,  the  vertical  degrees  of  freedom  can  be  computed  in  terms  of  the 
lateral  ones.  For  each  load  increment  the  kinematic  relations,  (6.12)  and  (6.13),  are  used 
as  predictors  to  remove  the  vertical  degrees  of  freedom  from  global  matrix  equation  of 
equilibrium.  They  are  used  again  as  correctors,  after  the  solution  of  the  reduced  system 
of  equations  has  been  achieved.  During  prediction  phase  x3,  X3,  Ax3,  Ax'3  are  computed 
using  (6.12)  and  (6.13).  In  addition,  during  the  prediction  phase  equation  (3.51)  are  used 
to  compute  xi,  xf,  Axi,  Ax'1;  x2,  x2,  Ax2,  Ax2.  During  the  correction  phase  x3,X3,  Ax3; 
Ax3  are  recomputed  using  (6.12)-(6.13)  and  all  deformation  dependent  matrices,  the 
equivalent  nodal  forces,  the  boundary  conditions  and  the  lengths  of  elements  are  cor- 
rected. Further  xi,  xl5  Axi,  Ax',  x2,  x2,  Ax2,  Ax2  are  recomputed  using  (3.51).  This 
predictor-corrector  scheme  is  repeated  until  convergence  is  achieved  for  each  load  incre- 
ment. The  algorithm  is  summarized  in  Figure  6.3. 
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_/v 


Assume  Riser  Initial  Position 
at  t = 0 


Apply  Initial  Velocity  V 
and  Calculate  Load  Vector  q(  V) 


Calculate  [M],  [C],  [K],  {f} 
and  Element  length 


Calculate  Axi,  Ax\,  Ax2,  Ax'2 
From  Eq.(3.51) 

l 


Calculated,  *3,  Ad,  Ax3 

From  Extensibility  Condition 


* t + At 


Use  d,  d>  d,  x2,  Axi,  Ad,  Ad 

From  Time  Step  t 


Calculate  Aq 


Calculate  F using  Local  Equilibrium 

— \l/ 

Calculate  Te  = F • 1 and  st  = Te/EA 


Calculate  [M],  [C],  [K],  {f} 
and  Element  length 

k 

Calculate  Ad,  Ax^,  Ax2,  Ax'2 
From  Eq.(3.51) 

k 

Correct  x3,x'3,  Ax 3,  Ax3 
using  Extensibility  Condition 


YES 


Convergence  Check 


NO 


Update  the  values  of 

d , d , d,  d 
7K 


Figure  6.3:  Solution  Algorithm  Employing  Predictor-Corrector  Scheme 
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6.3  Vortex-Induced  Vibration  Due  to  Current 
In  section  3.5,  the  matrix  equation  of  equilibrium  for  current-vortex  model  (Eqs. 
(3.6 1)-(3.63))  is  derived.  As  discussed  in  the  section,  Eqs.  (3.61)  and  (3.62)  represent  the 
self-excited  oscillator  equations,  while  (3.63)  presents  inline  vibration  due  to  current.  In 
those  equations,  the  system  matrix  My  and  My  represent  the  system  mass,  Cy  and  Cy 
coriolis  force  terms,  C/y  and  Cy  fluid  damping  terms,  Ty  and  Ty  stiffness  terms  due  to 

A ^ 

effective  tension,  Cey  and  Ceij  stiffness  terms  due  to  centrifugal  force,  Ky  and  Ky  stif- 
fness terms  due  to  bending  of  pipe,  Fy  and  Fy  self-excited  forcing  terms  due  to  the  vor- 
tex shedding  generated  by  inline  current,  and  My  and  Ky  represent  mass  and  stiffness 
term  of  fluid  oscillator  respectively.  And  also,  the  system  vector  f j represents  the  nonlin- 
ear drag  damping  due  to  current  and/ j represents  the  nonlinear  forcing  term  in  fluid  os- 
cillator equation.  The  matrices  representing  coriolis  force  term,  which  cause  the 
distortion  of  the  response  of  a riser,  are  skew  symmetric,  while  the  others  are  symmetric. 

From  the  inspection  of  coefficient  matrices  and  vectors  in  (3.64),  it  may  be  recog- 
nized that  the  drag  force  in  inline  vibration  equation  and  the  fluid  force  in  fluid  oscillator 
equation  have  nonlinear  term.  In  other  words,  there  is  no  highly  nonlinear  terms  in  the 
equation  to  be  solved,  so  that  simple  iteration  can  be  applied  for  the  convergence  of  solu- 
tion instead  of  predictor-corrector  scheme.  Thus,  the  application  of  Newmark  method 
combined  with  the  simple  iteration  is  enough  to  solve  the  system  described  in  section  3.5. 
Using  the  method  mentioned  above,  a finite  element  code,  which  is  called  CODEV,  is 
developed.  Overall  structure  of  the  code  is  the  same  as  that  of  CODE  given  in  Fig.  4.2 
but  some  of  the  routines  should  be  modified  with  slightly  different  arrays  and  variables 
because  of  the  different  coefficient  matrices  and  vectors  of  system  as  shown  in 
(3 .61  )-(3 . 64) . For  instance,  in  routine  ELEMM,  ELEMV,  GETCON,  new  array  and  vari- 
ables such  as  a hidden  flow  variable  w should  be  defined  throughout  CODEV.  And  also, 
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variable  CONS(J,I)  in  RCON,  one  of  input  routine,  need  more  constants  as  written 
below. 

CONS(J,I)  = input  constants  and  velocities;  I = element  number 

J = 1 ; the  velocity  of  internal  fluid 
2;  current  velocity  at  top 
3;  current  velocity  at  bottom 
4;  Cd 
5,  #4 
6;  p(x3) 

7;  a\ 

8;  a\ 

9;  (o  s 
10;  a'2 

11;  Ca 

J — 5 - 10;  See  section  2.5.4. 

Further,  routine  NEWMRK  should  be  coded  to  perform  time  step  integration  and 
simple  iteration  as  discussed  above.  Figure  6.4  shows  the  flow  charts  of  components  of 
subroutine  NEWMRK  which  completes  Fig.  4.2.  Routine  NEWMRK  performs  Newmark 
time  step  integration  and  calls  other  subroutines  within  the  routine.  Routine  FORMGV, 


Figure  6.4:  Flow  Charts  of  Components  of  Subroutine  NEWMRK  in  CODEV 
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APLBCV,  and  SOLVE  just  belong  to  the  iteration  phase,  while  the  entire  routine  are  op- 
erated for  every  time  step.  Routine  STRESS  calculates  the  internal  forces  and  moments, 
and  RARANG  rearranges  response  vectors  to  include  or  exclude  the  fixed  values  at 
boundary. 

For  the  verification  of  CODEV,  the  numerical  results  are  inspected  through  the 
physical  interpretation  with  the  generally  expected  trends.  The  same  material  properties 
and  constants  of  riser,  which  have  been  already  applied  in  section  5.3.2,  are  selected 
again  in  this  section.  The  first  concern  is  that  the  deflected  shape  of  riser  should  be  sym- 
metrical when  the  riser  system,  with  the  exclusion  of  internal  tension  and  internal  flow,  is 
subject  to  a uniform  current  in  the  water  column.  Table  6.1  shows  the  responses  of  riser 
at  certain  time  in  a uniform  current.  Nodal  displacements  and  rotations  are  shown  to  be 
exactly  symmetrical  or  skew-symmetrical  with  respect  to  middle  node  26  (typed  " * " in 
front  of  node  number),  that  is,  the  deflected  shapes  are  exactly  symmetrical  when  there 
are  no  internal  tension  and  flow.  Next,  including  internal  tension  and  internal  flow,  the 
response  of  riser  is  examined.  Figure  6.5  presents  maximum  transverse  displacement  ver- 
sus current  velocity.  Since  the  transverse  vibrations  due  to  current  are  controlled  by  the 
oscillating  vortex  generated  around  the  riser,  they  must  have  resonant  bands  near  the 
critical  current  velocity  in  which  the  vortex  shedding  frequency  accords  with  the  system 
frequency.  For  instance,  the  1st  system  frequency  is  already  obtained  from  Table  5.2. 
The  vortex  shedding  frequency  can  be  also  calculated  from  Eq.  (2.72).  Taking  Strouhal 
number  S as  0.20,  the  critical  velocity  can  be  obtained  from  the  following  equation: 

to,  = InSuc/Dg  = oo  i 
to.?  = 271(0.20  )wc/{(2.0)(0.3048)}  = 0.81701 


uc  - 0.40  m/sec 


(6.14) 
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Table  6.1:  Nodal  Displacements  and  Rotations  Due  to  Uniform  Current  in  Transverse 
and  Inline  Direction  ( Time  = 5 sec,  Top  Tension  = 0,  w = 0,  Internal  Velocity  = 0 ) 


Node 

# 

Transverse  Response 

Inline  Response 

Displacement  (m) 

Rotation  (rad) 

Displacement  (m) 

Rotation  (rad) 

1 

.00000D+00 

.83193D-05 

■00000D+00 

.78644D-02 

2 

.66594D-01 

.84558D-05 

.23945D-01 

.78390D-02 

3 

.13422D+00 

.84810D-05 

.47737D-01 

.77651D-02 

4 

.20041D+00 

.81095D-05 

.71236D-01 

.76467D-02 

5 

.26290D+00 

.7591  ID-05 

.94310D-01 

.74872D-02 

6 

.32099D+00 

.70087D-05 

.1 1684D+00 

.72902D-02 

7 

.37484D+00 

.65810D-05 

.13872D+00 

.70590D-02 

8 

.42619D+00 

.63184D-05 

.15984D+00 

.67973D-02 

9 

.47514D+00 

.60141D-05 

.18013D+00 

.65091D-02 

10 

.52311D+00 

.61228D-05 

.19950D+00 

.61985D-02 

18 

.77795D+00 

.29004D-05 

.31556D+00 

.32081D-02 

19 

.80231D+00 

.32003D-05 

.32472D+00 

.28022D-02 

20 

.82770D+00 

.3121 2D-05 

.33264D+00 

.23961D-02 

21 

.85125D+00 

.27309D-05 

.33933D+00 

.1991  ID-02 

22 

.87054D+00 

.21307D-05 

.34478D+00 

.15881D-02 

23 

.88521  D+00 

.15961D-05 

.34901  D+00 

.11874D-02 

24 

.89656D+00 

.12147D-05 

.35203D+00 

.78929D-03 

25 

.90425D+00 

.69518D-06 

.35383D+00 

.39377D-03 

*26 

.90711  D+00 

.14229D-16 

.35443D+00 

.1 1 152D-14 

27 

.90425D+00 

-.69518D-06 

.35383D+00 

-.39377D-03 

28 

.89656D+00 

-.12147D-05 

.35203D+00 

-.78929D-03 

29 

.88521D+00 

-.15961D-05 

.34901  D+00 

-.1 1874D-02 

30 

.87054D+00 

-.21307D-05 

.34478D+00 

-.15881D-02 

31 

.85125D+00 

-.27309D-05 

.33933D+00 

-.1991  ID-02 

32 

.82770D+00 

-.31212D-05 

.33264 D+00 

-.23961D-02 

33 

.80231  D+00 

-.32003D-05 

.32472D+00 

-.28022D-02 

34 

.77795D+00 

-.29004D-05 

.31556D+00 

-.32081D-02 

42 

.523  il  D+00 

-.61228D-05 

. 1 9950D+00 

-.61985D-02 

43 

.47514D+00 

-.60141D-05 

.18013D+00 

-.65091D-02 

44 

.42619D+00 

-.63184D-05 

.15984D+00 

-.67973D-02 

45 

.37484D+00 

-.65810D-05 

.13872D+00 

-.70590D-02 

46 

.32099D+00 

-.70087D-05 

.1 1684D+00 

-.72902D-02 

47 

.26290D+00 

-.7591  ID-05 

.94310D-01 

-.74872D-02 

48 

.20041  D+00 

-.81095D-05 

.71236D-01 

-.76467D-02 

49 

.13422D+00 

-.84810D-05 

.47737D-01 

-.77651D-02 

50 

.66594D-01 

-.84558D-05 

.23945D-01 

-.78390D-02 

51 

.00000D+00 

-.83193D-05 

.00000D+00 
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Figure  6.5:  Maximum  Vortex-Induced  Displacement  in  Transverse  Direction  According 
to  the  Variation  of  Uniform  Inline  Current.  The  System  Includes  Internal  Tension  and 
Internal  Flow.  ( Internal  Velocity  = 6.0  m/sec  ) 
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This  critical  velocity  is  consistent  with  the  velocity  in  which  the  first  peak  takes  place  as 
shown  in  Fig.  6.5.  In  the  same  manner,  the  critical  velocities  for  second  and  third  peak  in 
Fig.  6.5  are  compared  with  the  values  calculated  using  Eq.  (6.14)  and  found  to  be  consis- 
tent with  each  other.  While  the  loading  due  to  vortex  shedding  in  the  transverse  direction 
is  oscillating  well  to  induce  harmonic  riser  motion,  the  constant  load  due  to  current  force 
in  inline  direction  causes  a steady  static  response.  This  inline  displacement  simply  be- 
comes larger  when  the  current  velocity  increases.  Figure  6.6  presents  maximum  inline 
displacement  due  to  uniform  current  which  shows  the  trend  as  stated  above.  Figure  6.7 
plots  the  time  history  of  the  transverse  displacement  at  the  middle  point  of  the  riser.  In 
this  figure,  three  current  velocities,  over,  at  and  below  first  critical  velocity  are  chosen  to 
demonstrate  the  magnification  in  amplitude  of  oscillating  displacement  resulting  from  the 
current  velocity  that  generates  resonant  vorticies.  The  middle  figure  corresponds  to  the 
case  of  resonance.  As  discussed  already,  all  three  figures  represent  harmonically  oscillat- 
ing ones  due  to  the  vortex  shedding.  These  harmonic  responses  due  to  vortex  shedding 
are  calculated  by  forcing  riser  with  an  imposed  initial  force.  The  vortex  driven  oscillation 
usually  reaches  to  a steady  state  in  about  10-15  cycles  for  top-hinged  case.  Also,  the 
slight  modulation  in  amplitude  can  be  detected  from  the  figure.  This  modulation  results 
from  a third  power  term  in  the  hidden  flow  variable  w on  the  right  side  of  fluid  oscillator 
equation  (2.72).  This  can  be  proved  from  the  fact  that  this  modulation  would  disappear  if 
the  third  power  term  is  excluded  from  the  fluid  oscillator  equation.  Figure  6.8  plots  the 
time  history  of  inline  displacement  at  the  middle  of  riser.  As  expected,  the  displacement 
in  the  current  direction  approaches  to  a steady  static  state  after  a few  cyclic  fluctuations. 
Figure  6.9  plots  the  time-dependent  trajectories  of  displacement  vector  at  the  middle  of 
the  riser.  The  first  trajectory  represents  the  trajectory  for  a current  velocity  below  first 
critical  velocity,  the  second  is  for  the  first  critical  velocity,  and  the  third  is  for  the  case 
above  the  first  critical  velocity.  As  discussed  in  Fig.  6.5  through  6.8,  similar  differences 
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Figure  6.6:  Maximum  Inline  Displacement  in  Current  Direction  According  to  the  Varia- 
tion of  Uniform  Inline  Current.  The  System  Includes  Internal  Tension  and  Internal  Flow. 
( Internal  Flow  Velocity  = 6.0  m/sec  ) 
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Current  Velocity  = 0.30  m/sec 


Current  Velocity  = 0.46  m/sec 


Figure  6.7:  Transverse  Displacement  Due  to  Inline  Current  at  the  Middle  Node  of  Riser 
According  to  Change  of  Time.  The  System  Includes  Internal  Tension  and  Internal  Flow 
of  6.0  m/sec. 
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Figure  6.8:  Inline  Displacement  Due  to  Inline  Current  at  the  Middle  Node  of  Riser  Ac- 
cording to  Change  of  Time.  The  System  Includes  Internal  Tension  and  Internal  Flow  of 
6.0  m/sec. 
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Current  Velocity  = 0.30  m/sec 


Current  Velocity  = 0.40  m/sec 


Current  Velocity  = 0.46  m/sec 


Figure  6.9:  Time-Dependent  Trajectories  of  Displacement  Vector  Tip  at  the  Middle  Node 
of  Riser.  The  System  Includes  Internal  Tension  and  Internal  Flow  of  6.0  m/sec. 
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in  the  displacement  trend  between  transverse  direction  and  inline  direction  are  also  found 
in  Fig.  6.9.  It  can  also  be  recognized  that  the  trajectory  shape  changes  from  almost  closed 
oval  to  irregular  zigzag.  This  phenomenon  results  from  the  shortening  of  the  response 
period  due  to  the  increase  of  vortex  shedding  frequency  as  shown  in  Fig.  6.7. 

Finally,  the  deflected  shapes  of  the  riser  system  including  and  excluding  internal 
flow  and  tension  are  plotted  in  Figure  6.10.  Four  different  combinations  of  plots  are  giv- 
en in  transverse  and  current  direction.  The  deflected  shapes  are  drawn  at  time  =100  sec. 
the  coriolis  force  effect  is  seen  to  be  negligible  on  the  inline  displacement.  This  is  due  to 
the  fact  that  the  response  of  the  riser  system  in  current  direction  approaches  to  static 
steady  state  after  a few  oscillations  as  shown  in  Fig.  6.8.  In  other  words,  100  sec  is 
enough  time  for  the  system  to  reach  a static  state  such  that  time-dependent  coriolis  force 
effect  disappears.  Thus,  for  the  case  with  internal  flow  inside  the  riser  system  but  with  no 
internal  tension,  the  curvature  of  the  deflected  shape  in  the  current  direction  simply  in- 
creases because  the  centrifugal  force  due  to  internal  flow  leads  to  a decrease  on  the  stif- 
fness. On  the  other  hand,  the  riser  displacement  in  inline  direction  is  drastically  reduced 
if  internal  tension  is  applied  to  the  system.  It  is  natural  because  the  stiffness  term  in  Eq. 
(3.63)  has  to  be  increased  due  to  internal  tension.  However,  the  fact  stated  above  is  not 
always  true  for  the  response  in  the  transverse  direction  because  the  riser  motion  is  always 
oscillating  regardless  of  the  passage  of  time.  The  coriolis  force  and  the  centrifugal  force 
both  appear  to  have  an  effect  on  the  vortex-induced  vibration  of  the  riser  system  all  the 
time  whenever  the  internal  fluid  flows  inside  the  riser.  More  explanations  are  given  in 
chapter  7.  The  deflected  shapes  for  the  case  of  no  internal  tension  are  very  different  from 
those  with  internal  tension  because  the  existence  of  the  internal  tension  changes  the  sys- 
tem frequency.  Therefore,  the  riser  system  seems  to  behave  as  expected  from  physical 
point  of  view. 
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Figure  6.10:  Deflected  Shape  of  Riser  System  in  Transverse  and  Current  Direction  Ac- 
cording to  the  Change  of  System.  Four  Different  Combinations  Are  Considered  As  De- 
scribed in  Figure.  (Current  Velocity  - 0.40  m/sec,  Internal  Flow  Velocity  = 6.0  m/sec) 
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6.4  Vibration  Due  to  Current  or  Wave  Loading 
The  matrix  equilibrium  equation  (3.63),  which  represents  inline  vibration  due  to 
current,  can  be  utilized  to  obtain  the  responses  of  riser  under  simple  current  or  wave 
loading  because  the  inline  vibration  is  not  coupled  with  the  transverse  vibration  which  is 
self-excited.  In  this  case,  the  CODEV  is  applied  by  condensing  the  motion  in  the  trans- 
verse direction  and  some  of  routines  are  slightly  modified  in  case  of  wave  loading.  The 
unidirectional  wave  induced  particle  velocity  and  acceleration  are  obtained  in  subroutine 
RWAV  based  on  linear  wave  theory.  Furthermore  ELEMV,  which  calls  RWAVE,  should 
also  be  modified  to  include  the  forcing  term  due  to  wave  and  variable  CONS(J,I)  in  rou- 
tine RCON  should  be  changed  into  the  different  type  of  constants  as  written  below. 


CONS(J,I)  = input  constants  and  velocities;  I = element  number 

J — 1 ; the  velocity  of  internal  2;  current  velocity 
3;  current  direction  to  x-axis 
4;  Cd  5;  Ca 

6;  wave  height  7;  wave  number 

8;  water  depth 
9;  wave  direction  to  x-axis 
10;  gravity 

As  described  above,  CODEV  is  modified  and  given  as  CODEL. 

The  application  of  CODEL  for  the  dynamic  analysis  of  the  riser  system  subject  to 
simple  current  or  wave  loading  does  not  always  give  good  result  because  it  is  pro- 
grammed on  the  basis  of  linear  model  developed  in  section  3.5.  Since  the  geometric  non- 

linearities  due  to  large  structural  deflections  become  significant  particularly  for  long 
risers  over  1000  m,  the  predictor-corrector  scheme  described  in  section  6.2  should  be  im- 
plemented to  solve  the  incremental  form  of  the  matrix  equilibrium  equation.  The  algo- 
rithm helps  to  avoid  the  possibility  of  divergence  due  to  iterations  and  to  improve 
convergence  rate.  CODE  is  written  for  this  purpose  and  NEWMRK,  one  of  the  routines 
in  CODE,  accomplishes  the  solution  algorithm  employing  the  predictor-corrector  scheme 
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as  depicted  in  Fig.  6.3.  Figure  6.11  depicts  the  flow  chart  of  components  of  routine 
NEWMRK.  Routine  FORMGM  and  APLBCM  is  included  in  NEWMRK  because  the 
system  mass,  damping,  and  stiffness  matrices  are  dependent  on  the  riser  deformations 
which  changes  with  time  and  for  every  iteration  step.  Routine  EXTEN  calculates  the  ver- 
tical degrees  of  freedom  using  extensibility  conditions  (See  Eqs.  (6.12)  and  (6.13)).  In 
subroutine  NEWMRK,  routine  FORMGM  through  STRESS  are  included  within  the  it- 
eration loop  and  are  called  upon  for  every  time  step.  The  step  calculating  the  internal  ten- 
sion is  performed  within  the  routine  STRESS. 

Verification  on  CODEL  is  not  carried  out  as  it  is  only  slightly  modified  from  CO- 
DEV  keeping  the  same  basic  frame.  Instead,  the  CODE  version  is  tested  by  examining 
whether  the  results  are  reasonable.  For  the  purpose  of  investigating  the  validity  of 
CODE,  the  comparison  of  the  results  obtained  from  CODEL  with  one  from  CODE  is 
given.  It  is  well  known  from  previous  studies  that  the  nonlinear  effects  due  to  large  struc- 
tural deflections  and  nonlinear  drag  damping  become  significant  as  riser  length  increases 


* * * * # 


* : See  Figure  4.2 

# : See  Figure  6.4 


Figure  6.1 1:  Flow  Chart  of  Components  of  Subroutine  NEWMRK  in  CODE 
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or  internal  tension  decreases.  On  the  contrary,  the  nonlinear  effects  disappear  for  the  riser 
with  high  tension  and  short  length.  In  this  section,  CODE  and  CODEL  are  both  run  using 
the  same  riser  with  parameters  listed  in  Figure  5.2,  and  then  compared  with  each  other  to 
see  the  difference  between  the  linear  and  nonlinear  models  due  to  the  change  of  riser 
length.  The  flow  field  is  generated  by  a linear  wave  and  a uniform  current  at  right  angle 
to  each  other. 

Figure  6.12  and  6.13  plot  time  simulations  of  displacement  at  the  middle  of  the  ris- 
er in  current  and  wave  direction,  respectively.  Two  different  length  of  risers,  one  long 
and  one  short,  are  used  with  same  top  and  bottom  tension.  From  those  figures,  the  reduc- 
tion of  nonlinear  effect  on  the  short  one  can  be  seen.  Also,  from  those  time  simulation 
plots,  the  nonlinear  effect  in  the  long  riser  causes  change  of  the  system  frequency.  The 
displacement  envelopes  at  a specific  time  are  presented  in  Figure  6.14  and  6.15  for  cur- 
rent and  wave  direction,  respectively.  The  envelopes  in  current  direction  are  plotted  at 
time  = 200  sec  after  the  riser  reaches  the  static  state,  while  the  envelopes  in  wave  direc- 
tion are  plotted  at  an  arbitrary  time  when  the  displacement  at  the  middle  point  of  riser  is 
near  the  positive  peak.  The  difference  between  linear  and  nonlinear  models  grows  as  riser 
length  increases.  The  nonlinear  effect  appears  to  stiffen  the  riser  for  current  loading  but 
soften  it  for  wave  loading.  These  trends  are  also  found  from  previous  study  (Bemitsas 
and  Kokarakis,  1987).  Figure  6.16  shows  the  maximum  displacement  in  the  wave  direc- 
tion according  to  wave  frequency.  Like  a vortex-induced  vibration,  the  response  in  the 
wave  direction  also  shows  resonant  peaks  near  the  system  frequencies  but  the  magnitudes 
are  drastically  reduced  for  higher  frequency  wave.  This  phenomenon  seems  to  result 
from  the  fact  that  the  wave  velocity  distribution  in  water  depth  decays  faster  for  shorter 
frequency  waves.  Subsequently,  time  histories  of  the  displacement  at  the  middle  of  a 152 
m riser  in  wave  direction  are  presented  in  Fig.  6.17  for  three  wave  frequencies  near  first 
peak.  The  response  reach  a steady  oscillating  state  within  a few  cycle  and  shows  irregular 
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Figure  6.12:  Time  Simulations  of  the  Displacement  at  the  Middle  of  Riser  in  Current 
Direction.  Risers  with  Different  Length  Are  Subject  to  a Uniform  Current  (U  = 0.4 
m/sec)  and  a Linear  Wave  (Hw  = 6.1m  and  Tw  = 7.8  sec)  with  a Right  Angle  at  the 
Same  Time.  (Dashed  Line  = Nonlinear  Model,  Solid  Line  = Linear  Model) 
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Time  (sec) 


Figure  6.13:  Time  Simulations  of  the  Displacement  at  the  Middle  of  Riser  in  Wave 
Direction.  Risers  with  Different  Length  Are  Subject  to  a Uniform  Current  (U  = 0.4 
m/sec)  and  a Linear  Wave  (Hw  = 6.1  m and  Tw  = 7.8  sec)  with  a Right  Angle  at  the 
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Displacement  Envelope  (m)  Displacement  Envelope  (m) 

at  Time  = 200  sec  at  Time  = 200  sec 


Figure  6.14:  Displacement  Envelope  in  Current  Direction  under  Current  (U  = 0.4  m/sec) 
and  Wave  (Hw  = 6.1  m and  Tw  = 7.8  sec)  with  a Right  Angle.  (Dashed  Line  = Nonlinear 
Model,  Solid  Line  = Linear  Model) 
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Displacement  Envelope  (m)  Displacement  Envelope  (m) 

at  Time  = 190  sec  at  Time  =187  sec 


Figure  6.15:  Displacement  Envelope  in  Wave  Direction  under  Current  (U  = 0.4  m/sec) 
and  Wave  (Hw  = 6.1  m and  Tw  = 7.8  sec)  with  a Right  Angle.  (Dashed  Line  = Nonlinear 
Model,  Solid  Line  = Linear  Model) 
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Figure  6.16:  Maximum  Displacement  in  Wave  Direction  According  to  the  Change  of 
Wave  Frequency.  ( Hw  = 6.1  m ) 
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Figure  6.17:  Time  Simulations  of  the  Displacement  at  the  Middle  of  a 152  m Riser  in 
Wave  Direction.  (Hw  = 6.1  m) 
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modulation  in  amplitude  due  to  the  nonlinearities.  Time  trajectories  plotted  for  the  above 
three  cases  are  very  irregular  as  seen  in  Fig.  6.18;  this  irregularity  may  be  due  to  the 
combined  effect  of  current  and  wave  applied  on  the  same  time.  However,  if  the  current 
approaches  the  riser  at  a right  angle  to  wave  after  the  riser  response  due  to  wave  loading 
reached  a steady  oscillating  state,  the  irregular  shapes  should  change  to  that  shown  in 
Figure  6.9. 

Even  though  the  structural  damping  is  not  still  considered  in  the  model,  we  can 
easily  count  for  a Rayleigh  type  damping  such  as  given  in  Eq.  (6.15). 

[C]=A0[M]+Al  [K]  (6.15) 

where  Aa  and  Hi  are  Rayleigh  damping  coefficients. 

Rayleigh  damping  can  be  implemented  on  CODE  simply  by  adding  a few  line  in  the  end 
of  routine  ELEMM.  Considering  the  structural  damping,  time  history  and  trajectory  of 
displacement  at  the  middle  of  riser  in  both  direction  are  presented  in  Figure  6.19.  In  this 
figure,  the  riser  responds  to  current  at  a right  angle  to  wave  direction  after  a steady  oscil- 
lating state  in  wave  direction  has  been  reached  and  the  top  tension  ratio  (TR  = Top  Ten- 
sion / wL)  is  changed  from  2.6  to  1.2  for  the  reduction  of  stiffness  term.  Since  a high 
velocity  of  current  (U  = 1.5  m/sec)  and  a reduced  top  tension  (TR  = 1.2)  are  selected,  it 
can  be  seen  clearly  from  the  first  figure  of  Fig.  6.19  that  the  oscillation  of  riser  in  current 
direction  reaches  a static  state  within  2 to  3 cycles  with  the  longer  period.  In  addition,  it 
can  be  recognized  from  the  third  figure  of  Fig.  6.19  that  the  time  trajectory  shape  is 
transformed  into  more  regular  shape  because  of  a wave  with  relatively  low  period  (Tw  = 
4 sec)  in  this  application.  This  trend  of  time  trajectory  can  be  also  found  in  Mcnamara's 
study  (1988). 
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Figure  6.18:  Time  Trajectories  of  the  Displacement  at  the  Middle  of  a 152  m Riser  under 
Wave  (Hw  = 6.1  m)  and  Current  (U  = 0.4  m/sec)  with  a Right  Angle. 
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Figure  6.19:  Time  Simulations  in  Current  and  Wave  Direction  and  Time  Trajectory  at  the 
Middle  of  a 152  m Riser  under  Current  (U  = 1.5  m/sec)  and  a Linear  Wave  (Hw  = 6.1  m 
and  Tw  = 4 sec)  with  a Right  Angle.  Rayleigh  Damping  [C]  = 0.01  [M]  + 0.03  [K]  and 
Top  Tension  Ratio  TR  =1.2 
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6.5  Modification  of  System  Matrices  or  Vector  for  Top  Boundary  Condition 

As  mentioned  in  section  2.5.2,  the  top  boundary  condition  can  be  modeled  by  using 
translational  or  rotational  springs  to  provide  restoring  boundary  forces  or  moments.  In 
addition,  structures  like  a platform  or  an  anchor  can  be  modeled  as  if  mass  and  damper 
are  attached  at  the  top  of  riser.  These  concentrated  masses,  dampers,  and  springs  are  sim- 
ply added  to  the  corresponding  system  mass,  damping,  and  stiffness  matrix.  However, 
three  directional  components  of  top  tension  should  be  added  to  the  corresponding  system 
vector.  As  shown  in  Eq.  2.53  through  2.55,  these  direction  cosines  have  deformation- 
dependent  features.  Thus,  the  direction  cosines  should  be  included  within  the  iteration 
routine.  The  deformation  dependency  of  the  direction  cosines  lost  in  linear  model  and 
therefore,  the  vertical  top  tension  should  be  added  in  the  corresponding  degree  of  global 
stiffness  matrix.  By  the  principle  introduced  above,  the  modification  of  system  matrices 
is  performed  to  include  the  attached  properties  at  the  top  boundary  in  subroutine 
MODGMB. 

The  numerical  application  has  been  made  using  the  riser  with  a massive  structure 
and  mooring  system  attached  at  the  top.  The  mass,  damping,  and  spring  constants  are  ar- 
bitrarily selected  only  for  the  illustrative  purpose.  The  arbitrary  selection,  however,  does 
not  affect  the  nature  of  results  obtained.  The  uniform  current  is  applied  over  the  entire 
riser  at  a right  angle  to  wave  direction  after  the  wave-induced  vibration  reaches  a steady 
state.  The  time  simulations  of  tip  displacement  in  wave  and  current  direction  are  plotted 
in  the  first  and  second  figure  of  Fig.  6.20,  respectively  and  the  trace  of  riser  tip  on  a hori- 
zontal plane  is  plotted  in  the  third.  Since  the  strong  mass  and  damping  are  attached  at  the 
top  of  riser,  it  takes  long  for  the  riser  response  in  both  directions  to  reach  a steady  state. 
In  addition,  the  nonlinear  drift  in  the  wave  direction  can  be  found  even  though  a linear 
wave  is  applied  along  the  riser.  Figure  6.21  shows  the  displacement  envelopes  in  the  cur- 
rent direction.  It  can  be  recognized  from  this  figure  that  the  riser  response  reaches  a static 
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Displacement  in  Wave  Direction  (m) 


Figure  6.20:  Time  Simulations  in  Current  and  Wave  Direction  and  Time  Trajectory  at  the 
Top  of  a 152  m Riser  under  Current  (U  = 1.5  m/sec)  and  a Linear  Wave  (Hw  = 10  m and 
Tw  = 16  sec)  with  a Right  Angle.  Rayleigh  Damping  [C]  = 0.01  [M]  + 0.03  [K]  and  Top 
Tension  Ratio  TR  =1.2 
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Displacement  Envelope  (m) 


Figure  6.21:  The  Displacement  Envelopes  in  Current  Direction  at  a Specific  Time  for  a 
152  m Riser  under  Current  (U  = 1.5  m/sec)  and  a Linear  Wave  (Hw  = 10  m and  Tw  = 16 
sec)  with  a Right  Angle.  Rayleigh  Damping  [C]  = 0.01  [M]  + 0.03  [K]  and  Top  Tension 
Ratio  TR  =1.2 
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state.  However,  it  can  be  seen  from  Fig.  6.22  that  the  motion  of  top  does  not  reach  steady 
state  because  wave  period  is  relatively  short  to  the  period  of  the  riser  system.  In  addition, 
the  nonlinear  drift  in  the  wave  direction  can  be  also  found  even  though  a linear  wave  pass 
through  the  riser. 
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Displacement  Envelope  (m) 


Figure  6.22:  The  Displacement  Envelopes  in  Wave  Direction  at  a Specific  Time  for  a 
152  m Riser  under  Current  (U  = 1.5  m/sec)  and  a Linear  Wave  (Hw  = 10  m and  Tw  = 16 
sec)  with  a Right  Angle.  Rayleigh  Damping  [C]  = 0.01  [M]  + 0.03  [K]  and  Top  Tension 
Ratio  TR  =1.2 


CHAPTER  7 

THE  INVESTIGATION  OF  THE  EFFECT  OF  INTERNAL  FLOW 
ON  RISER  DYNAMICS 


In  the  previous  chapter,  the  verification  on  all  the  programs  developed  in  this  paper 
has  been  given.  In  this  chapter  the  effects  of  internal  flow  on  riser  dynamics  such  as  sys- 
tem frequencies,  vortex-induced  vibration,  and  vibrations  due  to  current  or  wave  loading 
are  investigated  by  applying  the  numerical  model.  The  investigations  are  carried  out  by 
inspecting  how  internal  flow  affects  the  kinematics  and  the  dynamics  of  a riser,  in  par- 
ticular, the  effects  of  the  forcing  terms  due  to  the  internal  flow,  namely,  the  centrifugal 
and  the  coriolis  forces.  In  addition,  the  effect  due  to  the  inclusion  of  the  geometric  and 
fluidic  nonlinearities  are  discussed. 

7.1  The  Effect  on  System  Frequency 

For  the  deep  sea  application,  the  riser  system  is,  in  general,  designed  to  be  nearly 
neutrally  buoyant.  Under  such  condition,  the  effective  weight,  the  difference  between 
riser  weight  and  local  net  buoyancy,  which  results  from  the  combined  fluid  action  of  ex- 
ternal and  internal  hydrostatic  pressure,  is  equal  to  zero,  that  is,  the  balance  between  riser 
weight  and  local  net  buoyancy  is  locally  maintained  for  the  entire  length  of  the  riser. 
The  basic  parameters  governing  the  property  of  the  riser  are  the  same  as  used  in  chapter 
5.  A change  is  made  on  the  magnitude  of  top  tension  and  riser  length.  The  system  fre- 
quencies were  computed  for  different  water  depths  with  six  different  internal  flow  condi- 
tions and  three  different  top  tensions.  Figure  7.1  shows  the  relationship  between  water 
depth  and  system  frequency  of  a neutrally  buoyant  riser  system,  i.e.,  w = 0.  The  ordinate 
is  the  ratio  of  the  1st  mode  frequency  to  the  natural  frequency  of  the  system  with  null  in- 
ternal flow.  As  can  be  seen  in  the  figure,  the  system  frequency  is  drastically  affected  by 
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Figure  7.1:  Internal  Flowing  Fluid  Effect  on  1st  Mode  System  Frequency  of  a Neutrally 
Buoyant  Riser  System.  The  internal  flow  velocity  changes  from  1.5  m/sec  to  12.0  m/sec 
and  mean  tension  (same  as  top  tension)  does  from  null  to  50  kN. 
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the  internal  flow  for  the  null  tension  case.  Finally,  riser  buckles  when  the  flow  velocity 
approaches  a critical  velocity.  The  increase  of  top  tension,  on  the  other  hand,  counters  the 
internal  flow  effect  as  seen  on  the  second  and  third  figures.  With  adequate  top  tension, 
the  system  frequency  approaches  constant  as  water  depth  increases.  This  is  because  the 
influence  of  the  stiffness  term  becomes  negligible.  Figure  7.2  shows  the  separate  effects 
of  the  two  flow-induced  force  components  centrifugal  and  coriolis  forces  on  a riser  sys- 
tem. The  case  is  one  of  the  cases  in  Figure  7.1,  with  a mean  tension  of  50  kN  and  an  in- 
ternal flow  velocity  of  12  m/sec.  From  Figure  7.2,  it  can  be  seen  that  the  main  effect  on 
system  frequency  is  due  to  centrifugal  force.  However,  even  though  for  cases  shown  the 
coriolis  force  seems  to  have  a minor  influence  on  system  frequency,  it  remains  unan- 
swered on  its  effect  on  system  dynamics  such  as  displacement  and  stress. 

7.2  Comparison  between  the  Effects  of  Centrifugal  and  Coriolis  Forces 
As  discussed  earlier,  internal  flow,  as  it  travels  along  the  curved  path  in  the  riser, 
generates  centrifugal  and  coriolis  force.  These  dynamic  forces  exerted  against  the  riser, 
in  turn,  affect  the  dynamic  behavior  of  a riser.  The  centrifugal  force  reduces  the  stiffness 
of  a riser  whereas  the  coriolis  force  causes  the  dynamic  coupling  with  other  forces.  Their 
effects  on  riser  behavior  are  separately  examined  here.  Figure  7.3  and  7.4  present  the  dis- 
placement shapes  of  a riser  at  four  different  arbitrarily  selected  times  to  illustrate  how  the 
deflected  shape  changes  with  the  inclusion  of  the  centrifugal  force  only.  In  both  cases  the 
riser  is  hinged  at  the  top.  Figure  7.3  shows  the  case  under  a 5 sec  excitation  wave  and 
Figure  7.4  shows  the  case  under  a 10  sec  excitation  wave.  In  comparison  with  the  no  in- 
ternal flow  situation,  the  magnitude  of  the  displacement  envelopes  are  changed  but  the 
shapes  remain  undistorted.  In  other  words,  the  ratio  of  magnification  remains  constant 
along  the  riser.  Also,  it  should  be  noted  here,  since  centrifugal  force  reduces  the  stif- 
fness of  the  riser,  the  amplitude  of  the  riser  displacement  also  increases  with  increasing 
centrifugal  force.  However,  a phase  difference  is  introduced  into  the  system.  Therefore, 
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Figure  7.2:  Comparison  Between  the  Effects  of  Each  Forcing  Terms  ( Centrifugal  and 
Coriolis  Forces)  Due  to  Internal  Flow  on  1 st  System  Frequency.  Mean  Tension  = 50  kN, 
Internal  Flow  Velocity  = 12  m/sec 
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Figure  7.3:  The  Effect  of  Centrifugal  Force  on  Displacement  Envelope  of  Riser  at  an  Ar- 
bitrary Time.  (TR  = 1.2,  w = 3.86  kN/m,  Hw  = 6.1  m,  Tw=  15  sec) 
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Figure  7.4:  The  Effect  of  Centrifugal  Force  on  Displacement  Envelope  of  Riser  at  an  Ar- 
bitrary Time.  (TR  = 1.2,  w = 3.86  kN/m,  Hw  = 6.1  m,  Tw  = 5 sec) 


127 


the  amplitude  amplification  will  not  be  manifested  at  every  instance. 

Unlike  the  effect  of  centrifugal  force,  the  coriolis  force  alters  the  displacement 
shape  as  shown  in  Fig.  7.5  and  7.6.  As  can  be  seen,  the  magnitudes  of  displacement 
change  at  upper  and  lower  part  of  riser  are  different.  This  can  be  understood  from  the 
fact  that  the  coriolis  force  term  is  composed  of  a mixed  derivative  with  respect  to  time 
and  space  or  the  coriolis  matrix  in  matrix  equilibrium  equation  is  skew  symmetric.  The 
coriolis  force  alternately  couples  with  inertia  and  stiffness  according  to  the  passage  of 
time.  In  addition,  it  can  be  recognized  from  the  comparison  of  figure  7.5  with  7.6  that  its 
effect  is  pronounced  for  shorter  period  wave  excitations.  It  results  from  the  fact  that  the 
coriolis  force  is  a function  of  rotational  velocity  which  is  faster  under  shorter  period 
waves. 

In  summary,  the  centrifugal  force,  which  depends  only  on  the  curvature  of  riser 
deflection,  does  not  alter  the  displacement  shape,  while  the  coriolis  force,  which  has  a 
mixed  derivative  of  time  and  space,  distorts  the  displacement  shape.  Further,  it  can  be 
inferred  from  the  figures  in  this  section  that  the  system  frequencies  are  affected  by  both 
forces. 

7.3  The  Effect  on  Vortex-Induced  Vibration 

As  described  in  the  previous  chapter,  Vortex-induced  vibration  due  to  inline  cur- 
rent has  been  formulated  by  modifying  Iwan-Blevin's  model  to  include  the  effect  of  in- 
ternal flow  and  a computer  program  CODEV  has  been  developed  to  solve  it.  The 
transverse  vibration  due  to  vortex  shedding  with  the  inclusion  of  internal  flow  is  investi- 
gated here,  while  the  inline  vibration  in  current  direction  will  be  discussed  in  section  7.5. 

The  time  simulation  of  the  displacement  at  the  middle  node  of  a riser  is  presented 
in  Figure  7.7.  The  presence  of  internal  flow  in  this  case,  which  has  a very  high  top  ten- 
sion, causes  a slight  increase  in  amplitude  and  a slight  shift  in  phase.  Figure  7.8  presents 
the  trajectory  of  displacement  at  this  middle  node  according  to  the  lapse  of  time.  The 
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Figure  7.5:  The  Effect  of  Coriolis  Force  on  Displacement  Envelope  of  Riser  at  an  Arbi- 
trary Time.  (TR  = 1.2,  w = 3.86  kN/m,  Hw  = 6.1  m,  Tw  = 15  sec) 
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Figure  7.6:  The  Effect  of  Coriolis  Force  on  Displacement  Envelope  of  Riser  at  an  Arbi- 
trary Time.  (TR=  1.2,  w = 3.86  kN/m,  Hw  = 6.1  m,  Tw  = 5 sec) 


Displacement  (m)  Displacement  (m) 


130 


V = 0.0  m/sec 


Figure  7.7:  The  Effect  of  Internal  Flow  on  Displacement  at  the  Middle  of  152  m Riser. 
(Current  Velocity  U = 0.32  m/sec,  Top  Tension  = 1550  kN,  The  Other  Properties  are  the 
Same  as  Example  of  Fig.  7.6) 
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Figure  7.8:  The  Effect  of  Internal  Flow  on  the  Trajectory  at  the  Middle  of  152  m Riser 
According  to  the  Lapse  of  Time.  (Current  Velocity  U = 0.32  m/sec,  Top  Tension  = 1550 
kN,  The  Other  Properties  are  the  Same  as  Example  of  Fig.  7.6) 
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trajectory  shape  is  modified  more  so  in  the  transverse  direction  than  the  inline  direction. 
This  is  partly  due  to  the  high  tension  applied  at  top  of  riser.  As  top  tension  eases,  the  ef- 
fect of  internal  flow  on  the  riser  responses,  particularly  on  inline  displacement,  should 
become  more  pronounced. 

The  optimal  riser  tension  can  be  determined  by  minimizing  the  combined  stress  of 
the  axial  stress  resulting  from  internal  tension  and  the  bending  stress  from  lateral  deflec- 
tion. Thus,  although  high  top  tension  tends  to  reduce  the  riser  displacement  and  the  asso- 
ciated bending  stress,  there  is  a limit  to  its  application  because  the  increase  in  tensile 
stress  may  cut  the  benefit  on  the  reduction  in  bending  stress.  Since  direct  stress  increases 
with  water  depth,  top  tension  has  to  be  eased  for  deep  water  application.  It  is  known  from 
previous  investigations  that  the  effect  of  internal  flow  on  riser  dynamics  is  not  significant 
under  high  top  tension.  However,  the  question  is  whether  the  internal  flow  effects  should 
be  considered  for  deep  water  case  for  the  reason  give  above. 

For  comparison  purpose,  the  same  riser  in  a 152  m of  water  depth  as  used  in  pre- 
vious examples  but  with  very  low  top  tension  is  used  here.  Computations  are  then  per- 
formed for  three  types  of  tension  (high,  middle,  and  low),  each  with  three  internal 
velocities  (null,  high,  and  very  high).  Figure  7.9  shows  the  effect  of  internal  flow  on 
maximum  displacement  in  the  transverse  direction,  i.e.,  the  direction  that  the  forces  due 
to  vortex  shedding  act.  Since  the  vortex  model  proposed  by  Iwan  and  Blevin  has  proven 
to  work  well  in  certain  range  of  Reynolds  number,  current  velocities  selected  in  the  pres- 
ent computations  are  also  determined  from  equation  (2.72)  within  this  Reynolds  number 
range.  As  shown  in  the  figure,  the  presence  of  internal  flow  causes  only  a slight  shift  of 
the  resonant  bands  as  well  as  a slight  increase  in  peak  values.  However,  for  non-peak 
positions,  the  effect  of  internal  flow  on  displacement  amplitude  is  much  larger.  Follow- 
ing the  curve  of  any  set  one  sees  that  internal  flow  causes  both  amplitude  amplification 
and  reduction  depending  upon  the  flow  velocity.  This  is  clearly  significant  in  riser  design 
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Figure  7.9:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  in  Vortex  Direction 
for  152  m Riser.  Riser  Properties  except  Top  Tension  are  the  Same  as  Example  Figure 
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in  determining  the  optimum  internal  flow  velocity.  In  other  words,  displacement  ampli- 
tudes are  either  augmented  or  reduced  depending  on  where  the  current  velocity  magni- 
tudes as  shown  in  this  figure.  This  situation  is  also  shown  in  Figures.  7.10  and  7.11. 
Also,  as  internal  current  velocity  increases,  the  modal  shape  will  also  shift  from  lower  to 
higher  modes.  For  instance,  for  the  case  presented  in  Figure  7.10,  the  riser  deflection 
shape  changes  from  1st  to  2nd  mode  when  U = 0.6  m/sec  and,  again,  changes  to  3rd 
mode  when  U =0.9  m/sec. 

In  summary,  the  effect  of  internal  flow  on  vortex-induced  vibration  is  not  negli- 
gible, in  particular,  when  the  top  tension  is  low  as  the  cases  found  in  deep  water. 

7.4  Consideration  of  Geometric  Nonlinearity 

The  nonlinearities  in  riser  dynamic  analysis  are  mainly  of  two  origins  that  from  the 
flow-induced  drag  force  and  from  geometric  nonlinearities  due  to  large  structural  deflec- 
tions. In  particular,  the  latter  becomes  significant  for  long  riser.  To  examine  the  effects  of 
geometric  nonlinearity  maximum  displacement  are  of  riser  with  internal  flow  can  be 
computed  by  the  program  for  different  riser  lengths  or  water  depths.  But,  there  is  no 
readily  available  information  on  selecting  riser  parameters  for  every  water  depth  and  on 
determining  top  tension  value  to  avoid  riser  failure  due  to  excessive  stress.  A conserva- 
tive approach  is  taken  here  to  let  the  top  tension  maintained  a constant  regardless  of  wa- 
ter depth  and  then,  let  the  effective  weight  be  adjusted  to  avoid  negative  tension  over  the 
riser  length  for  different  water  depths.  In  practice,  this  effective  weight  can  be  adjusted 
by  altering  the  buoyancy  along  the  riser. 

With  this  hypothetical  riser,  the  maximum  displacements  of  the  riser  are  computed 
for  water  depth  up  to  1500  m using  both  linear  and  nonlinear  models.  The  ratios  of  these 
maximum  displacements  to  that  of  the  case  without  internally  flow  are  plotted  in  Fig- 
ure 7.12.  As  shown  in  the  figure,  this  maximum  displacement  ratio  basically  increases 
with  depth  to  a certain  peak  point  beyond  which  it  decreases  with  depth  and  gradually 
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Figure  7.10:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  Envelope  in  Vortex 
Direction  for  152  m Riser  at  Some  Current  Velocity.  Top  Tension  = 1550  kN,  Bottom 
Tension  = 1200  kN.  Riser  Properties  except  Tension  are  the  Same  as  Example  of  Figure 
7.6.  Top  Condition  is  Semi-Restrained. 
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Figure  7.11:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  Envelope  in  Vortex 
Direction  for  152  m Riser  at  Some  Current  Velocity.  Top  Tension  = 350  kN,  Bottom 
Tension  = 0 kN.  Riser  Properties  except  Tension  are  the  Same  as  Example  of  Figure 
7.6.  Top  Condition  is  Semi-Restrained. 
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Figure  7.12:  Comparison  Between  Linear  and  Nonlinear  Model  about  the  Effect  of  Inter- 
nal Flow  According  to  the  Change  of  Water  Depth  and  Internal  Flow  Velocity.  Top  Ten- 
sion = 650  kN,  TR  = 1.1,  Effective  Weight  = Top  Tension  / (TR  x Riser  Length) 
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approaches  a constant  as  the  effect  of  the  riser  stiffness  also  diminishes  with  depth.  The 
effect  due  to  nonlinear  geometry  can  also  be  examined  from  this  Figure.  As  can  be  seen, 
in  shallow  water  the  effect  of  geometric  nonlinearity  is  actually  to  reduce  the  displace- 
ment ratio.  However,  as  the  water  becomes  deeper  than  certain  critical  value,  the  dis- 
placement ratio  becomes  larger  with  the  consideration  of  the  geometrical  nonlinearity. 
The  differences  from  the  linear  model  and  the  nonlinear  model  clearly  are  not  negligible 
and  they  become  larger  for  higher  internal  flow  velocity.  For  the  case  of  increased  top 
tension  such  as  shown  in  Figure  7.13,  the  displacement  ratio  is  also  found  to  follow  the 
same  trend  but  the  differences  between  the  linear  and  nonlinear  models  become  smaller. 

In  summary,  the  effect  of  geometrical  nonlinearity  becomes  increasingly  more  im- 
portant for  risers  in  deeper  water,  particularly  if  the  top  tension  is  low  and  the  internal 
flow  velocity  is  high.  Under  such  conditions,  linear  model  underestimates  the  displace- 
ment as  compared  to  nonlinear  model.  A linear  model  can  be  used  if  top  tension  is  high 
which  generally  restricts  its  application  to  shallow  water. 

7.5  Riser  Vibrations  Due  to  Wave  or  Current  Loading 

In  this  section  the  system  vibrations  due  to  current  or  wave  loading  of  a riser  with 
internal  flow  are  examined.  The  CODE  program  developed  in  this  investigation  is  used 
to  perform  the  analysis.  As  mentioned  in  earlier  the  CODE  program  includes  the  nonlin- 
ear effects  due  to  flow-induced  drag  force  and  large  structural  deformation,  both  of 
which  become  increasing  important  in  deeper  waters. 

First,  the  time  simulations  of  tip  displacements  are  calculated  for  the  condition  of 
a combined  uniform  current  and  a wave  system  at  right  angle  to  each  other.  The  results 
are  shown  in  Figures  7.14  to  7.16.  Figure  7.14  shows  the  current-direction  inline  dis- 
placement which  behaves  like  a damped  oscillations  that  eventually  will  approach  a 
steady  displacement.  Since  a substantial  mass  with  strong  damping  and  spring  character- 
istics is  attached  at  the  top  of  riser  to  model  a massive  anchored  platform,  the  effect  of 
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Figure  7.13:  Comparison  Between  Linear  and  Nonlinear  Model  about  the  Effect  of  Inter- 
nal Flow  According  to  the  Change  of  Water  Depth  and  Internal  Flow  Velocity.  Top  Ten- 
sion = 765  kN,  TR  = 1.1,  Effective  Weight  = Top  Tension  / (TR  x Riser  Length) 
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Figure  7.14:  The  Effect  of  Internal  Flow  on  the  Displacement  in  Current  Direction  at  the 
Tip  of  152  m Riser.  Second  Figure  Shows  the  Zoom  Picture  at  Certain  Time  Interval. 
(Current  Velocity  U = 0.40  m/sec,  Top  Tension  = 350  kN,  Effective  Weight  w = 2.3 
kN/m,  Semi-Restrained  Top  Condition) 
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Figure  7.15:  The  Effect  of  Internal  Flow  on  the  Displacement  in  Wave  Direction  at  the 
Tip  of  152  m Riser.  Second  Figure  Shows  the  Zoom  Picture  at  Certain  Time  Interval. 
(Wave  Period  Tw  = 20  sec,  Wave  Height  Hw  = 6.1  m,  Top  Tension  = 350  kN,  Effective 
Weight  w = 2.3  kN/m,  Semi-Restrained  Top  Condition) 
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Displacement  in  Current  Direction  (m) 


Figure  7.16:  The  Effect  of  Internal  Flow  on  the  Time-Dependent  Trajectory  at  the  Tip  of 
152  m Riser.  (Current  Velocity  U = 0.4  m/sec,  Wave  Period  Tw  = 20  sec,  Wave  Height 
Hw  = 6.1  m,  Top  Tension  = 350  kN,  Effective  Weight  w = 2.3  kN/m,  Semi-Restrained 
Top  Condition,  Current  and  Wave  Cross  with  an  Right  Angle) 
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internal  flow  is  clearly  overshadowed  in  this  case.  Figure  7.15  plots  the  time  displace- 
ment curve  in  the  wave  direction.  It  is  shown  the  effect  of  internal  flow  is  more  pro- 
nounced in  this  direction  owing  to  the  fact  that  the  coriolis  force  is  maintained  by  the 
wave  oscillation.  Figure  7.16  shows  the  time-dependent  trajectory  of  riser's  tip  displace- 
ment. As  expected,  the  coupled  response  in  both  directions  results  in  the  narrowing  pat- 
tern in  the  current  direction  but  maintains  a modulated  oscillation  in  the  wave  direction. 
The  presence  of  internal  flow  causes  a slight  decrease  in  current-direction  displacement 
but  a more  pronounced  increases  in  wave-direction  displacement.  Figure  7.17  presents 
the  effect  of  internal  flow  on  the  maximum  displacement  as  a function  of  current  strength 
for  different  riser  tensions  and  internal  velocities.  In  this  case,  the  centrifugal  force  domi- 
nates the  coriolis  effect  and  the  effect  of  internal  flow  becomes  larger  with  increasing 
current  velocity  as  a result  of  the  increased  riser  deformation.  This  effect  is  farther  mag- 
nified for  risers  with  low  tension  as  seen  in  this  figure.  In  other  words,  in  order  to  reduce 
the  centrifugal  force  due  to  internal  flow  one  must  increase  either  the  top  tension  or  the 
buoyancy.  These  effects  also  can  be  found  in  Figures  7.18  and  7.19  which  present,  re- 
spectively, the  cases  of  internal  flow  on  maximum  envelope  of  riser  under  high  and  low 
top  tensions.  Under  high  top  tension,  the  centrifugal  force  effect  is  relatively  insignifi- 
cant, while  under  low  tension  its  effect  becomes  more  important  as  current  velocity  in- 
creases. Figure  7.20  plots  the  maximum  displacement  in  the  wave  direction  as  a function 
of  wave  frequency  for  different  top  tensions  and  internal  flow  velocities.  Like  the  case  of 
vortex-induced  vibration,  displacement  due  to  wave  loading  also  has  resonance  peaks  of 
which  the  positions  are  lightly  shifted  due  to  the  presence  of  internal  flow.  The  first  reso- 
nance peak  apparently  dominates  the  rest  for  all  the  cases.  Here,  the  fundamental  natural 
frequency  of  the  low  tension  case  (about  0.4  rad/sec  in  Fig.  7.22)  is  approximately  half 
that  of  the  high  tension  case  (about  0.8  rad/sec  in  Fig.7.21).  In  these  first  mode  frequency 
bands,  coriolis  force  appears  to  be  more  important  than  the  centrifugal  force  in  internal 
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Figure  7.17:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  According  to  the 
Increase  of  Current  Velocity.  System  is  Same  as  Example  of  Fig.  7.16  except  Parameters 
Shown  in  this  Figure. 
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Figure  7.18:  The  Effect  of  Internal  Flow  on  Maximum  Envelopes  of  a 152  m Riser  Ac- 
cording to  the  Change  of  Current  Velocity  under  High  Top  Tension  of  1550  kN.  System 
is  Same  as  Example  of  Fig.  7.16  except  Parameters  as  Shown  in  This  Figure. 
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Figure  7.19:  The  Effect  of  Internal  Flow  on  Maximum  Envelopes  of  a 152  m Riser  Ac- 
cording to  the  Change  of  Current  Velocity  under  Low  Top  Tension  of  350  kN.  System  is 
Same  as  Example  of  Fig.  7.16  except  Parameters  as  Shown  in  This  Figure. 
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Figure  7.20:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  According  to  the 
Increase  of  Wave  Frequency.  System  is  Same  as  Example  of  Fig.  7.16  except  Parameters 
Shown  in  this  Figure. 
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Figure  7.21:  The  Effect  of  Internal  Flow  on  Maximum  Envelopes  of  a 152  m Riser  Ac- 
cording to  the  Change  of  Wave  Frequency  under  High  Top  Tension  of  1550  kN.  System 
is  Same  as  Example  of  Fig.  7.16  except  Parameters  as  Shown  in  This  Figure. 
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Figure  7.22:  The  Effect  of  Internal  Flow  on  Maximum  Envelopes  of  a 152  m Riser  Ac- 
cording to  the  Change  of  Wave  Frequency  under  Low  Top  Tension  of  350  kN.  System  is 
Same  as  Example  of  Fig.  7.16  except  Parameters  as  Shown  in  This  Figure. 
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flow  effects.  This  is  observed  from  the  reduced  peak  values  due  to  the  positive  effect  of 
coriolis  force.  Centrifugal  force  should  have  the  adverse  effect  of  increasing  the  displace- 
ment. The  same,  however,  can  not  be  said  for  higher  modes.  Figures  7.21  and  7.22  also 
shows  that  the  internal  flow  effect,  be  it  coriolis  or  centrifugal,  plays  a more  important 
role  under  low  top  tension  as  evidenced  from  comparing  the  changes  of  modal  envelopes 
due  to  changes  of  internal  velocities  in  these  two  cases. 

In  summary,  the  effect  of  internal  flow  on  tip  displacement  is  not  significant  if  the 
riser  is  restrained  at  the  top  by  a large  mass  with  strong  damping  and  stiff  spring.  The 
effect  of  internal  flow  on  the  displacement  of  a riser  in  a current  field  is  dominated  by 
centrifugal  force  because  the  riser  response  is  slowly  time  varying  but  could  have  large 
curvature.  On  the  other  hand,  for  riser  under  pure  wave  loading  the  internal  flow  effect  is 
dominated  by  coriolis  force  owing  to  the  oscillatory  nature  of  riser  response.  Also,  cen- 
trifugal force  always  increases  the  displacement  which  is  undesirable,  but  coriolis  force 

could  induce  reduced  displacement  at  certain  resonant  peaks. 

7.6  Effects  on  Internal  Stresses  and  Displacement  Fields 

In  this  section,  the  effects  of  internal  flow  on  riser  kinematics  and  dynamics  such 

as  displacement,  rotation,  shear  force,  and  bending  moment,  are  investigated.  The  re- 
sponses due  to  current  loading  and  vortex  shedding  are  addressed  specifically.  The  re- 
sponses due  to  wave  loading  are  of  same  nature  as  due  to  vortex  shedding. 

As  shown  in  Fig.  7.23,  the  effects  of  internal  flow  on  maximum  displacement  and 
rotation  become  more  important  when  current  velocity  and/or  internal  flow  velocity  in- 
crease, and  when  top  tension  decreases.  In  particular,  it  can  be  seen  from  Fig.  7.24  that 
the  maximum  effect  on  displacement  occurs  at  locations  where  the  slopes  are  zero  and 
the  maximum  effect  on  riser  rotation  occurs  at  the  bottom.  The  effects  on  bending  mo- 
ment are  the  same  as  on  displacement  but  the  effects  on  shear  force  are  not  as  evidenced 
in  Figures  7.25  and  7.26.  For  the  system  under  harmonic  loading  such  as  vortex  shedding 
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Figure  7.23:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  Fields  According  to 
the  Increase  of  Current  Velocity.  System  is  Same  as  Example  of  Fig.  7.16  except  Pa- 
rameters Shown  in  this  Figure. 
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Figure  7.24:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  and  Rotation  Enve- 
lope Due  to  Current  Loading  for  Three  Types  of  Top  Tension.  System  is  Same  as  Exam- 
ple of  Fig.  7.16  except  Parameters  Shown  in  this  Figure. 
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Figure  7.25:  The  Effect  of  Internal  Flow  on  Maximum  Internal  Stress  Fields  Due  to  Cur- 
rent According  to  the  Increase  of  Current  Velocity.  System  is  Same  as  Example  of  Fig. 
7.16  except  Parameters  Shown  in  this  Figure. 
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Figure  7.26:  The  Effect  of  Internal  Flow  on  Maximum  Shear  and  Bending  Moment  En 
velope  Due  to  Current  Loading  for  Three  Types  of  Top  Tension.  System  is  Same  as  Ex 
ample  of  Fig.  7.16  except  Parameters  Shown  in  this  Figure. 
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or  wave,  the  resonance  bands  of  both  displacement  and  rotation  are  affected  slightly  and 
the  changes  become  more  visible  under  low  tension  (See  Fig.  7.27).  Also,  it  can  seen 
from  Fig.  7.28  the  maximum  effect  on  displacement  is  at  locations  where  the  slope  is 
zero,  while  the  maximum  effect  on  rotation  occurs  not  only  at  bottom  but  also  at  inflec- 
tion points.  Shear  force  or  bending  moment  also  experience  changes  of  resonant  bands  in 
the  same  fashion  as  the  rotation  or  displacement  does  (See  Fig.  7.29).  Furthermore, 
maximum  effect  on  shear  force  occurs  at  inflection  points,  as  shown  in  Figure.  7.30. 
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Figure  7.27:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  Fields  Due  to  Vor- 
tex Shedding  According  to  the  Increase  of  Current  Velocity.  System  is  Same  as  Example 
of  Fig.  7.16  except  Parameters  Shown  in  this  Figure. 
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Figure  7.28:  The  Effect  of  Internal  Flow  on  Maximum  Displacement  and  Rotation  Enve 
lopes  Due  to  Vortex  Shedding  for  Three  Types  of  Top  Tensions.  System  is  Same  as  Ex 
ample  of  Fig.  7.16  except  Parameters  Shown  in  this  Figure. 
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Figure  7.29:  The  Effect  of  Internal  Flow  on  Maximum  Stress  Fields  Due  to  Vortex  Shed- 
ding According  to  the  Increase  of  Current  Velocity.  System  is  Same  as  Example  of  Fig. 
7.16  except  Parameters  Shown  in  this  Figure. 
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Figure  7.30:  The  Effect  of  Internal  Flow  on  Maximum  Stress  Envelopes  Due  to  Vortex 
Shedding  for  Three  Types  of  Top  Tensions.  System  is  Same  as  Example  of  Fig.  7.16  ex- 
cept Parameters  Shown  in  this  Figure. 


CHAPTER  8 

CONCLUSIONS  AND  FURTHER  STUDY 


8.1  Conclusions 

A mathematical  model  for  the  analysis  of  riser  systems  with  the  inclusion  of  inter- 
nal flow  is  developed  together  with  numerical  solutions  and  computer  programs.  Four 
different  programs  are  developed  in  this  study;  FRISER  program  computing  the  system 
frequencies  by  using  semi-analytical  method,  CODE  solving  a nonlinear  model  by  imple- 
menting predictor-corrector  scheme,  CODEV  computing  vortex-induced  vibrations,  and 
CODEL  obtaining  the  solutions  of  a linear  model.  From  the  results  of  sample  computa- 
tions, the  effects  of  internal  flow  on  riser  dynamics  are  examined.  The  following  conclu- 
sions are  drawn: 

1)  There  are  two  dynamic  forces  due  to  the  motion  of  internal  fluid,  that  is,  cen- 
trifugal and  coriolis  force.  The  centrifugal  force,  which  depends  only  on  the  curvature  of 
riser  deflection,  does  not  alter  displacement  shape  whereas  the  coriolis  force,  which  is  a 
term  with  mixed  derivative  of  time  and  space,  distorts  the  displacement  shape. 

2)  The  system  frequency  is  significantly  affected  by  the  internal  flow  for  the  null 
tension  case.  The  riser  will  buckle  finally  when  flow  velocity  reaches  a critical  velocity. 
The  increase  of  top  tension,  on  the  other  hand,  partially  counters  the  internal  flow  effects. 
Further,  the  main  effect  on  system  frequency  is  due  to  centrifugal  force. 

3)  The  effect  of  internal  flow  on  vortex-induced  vibration  is  not  negligible,  in  par- 
ticular, when  the  top  tension  is  low  as  for  the  cases  in  deep  water. 

4)  The  effect  of  geometrical  nonlinearity  becomes  increasingly  more  important  for 
risers  in  deeper  water,  particularly  if  the  top  tension  is  low  and  the  internal  flow  velocity 
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is  high.  Under  such  conditions,  the  linear  model  underestimates  the  displacement  as  com- 
pared to  the  nonlinear  model.  A linear  model  can  be  used  if  top  tension  is  high  which 
generally  restricts  its  application  to  shallow  water. 

5)  The  effect  of  internal  flow  on  tip  displacement  is  not  significant  if  the  riser  is 
restrained  at  the  top  by  a large  mass  with  strong  damping  and  a stiff  spring. 

6)  The  effect  of  internal  flow  on  the  displacement  of  a riser  in  a current  field  is 
dominated  by  centrifugal  force.  On  the  other  hand,  for  riser  under  pure  wave  loading  the 
internal  flow  effect  is  dominated  by  coriolis  force.  Also,  centrifugal  force  always  in- 
creases the  displacement  which  is  undesirable,  but  coriolis  force  could  increase  or  reduce 
displacement  depending  on  the  modes  of  resonant  peaks. 

7)  The  effects  of  internal  flow  on  maximum  displacement  and  rotation  become 
more  important  when  current  velocity  and/or  internal  flow  velocity  increase,  and  when 
top  tension  decreases.  In  particular,  the  maximum  effect  on  displacement  occurs  at  loca- 
tions where  the  slopes  are  zero  and  the  maximum  effect  on  riser  rotation  occurs  at  the 
bottom.  The  effects  on  bending  moment  are  the  same  as  on  displacement  but  the  effects 
on  shear  force  are  not  in  certain  rule. 

8)  For  the  system  under  harmonic  loadings  such  as  vortex  shedding  or  wave,  the 
resonance  bands  of  both  displacement  and  rotation  are  affected  slightly  and  the  changes 
become  more  visible  under  low  tension.  Also,  the  maximum  effect  on  displacement  is  at 
locations  where  the  slope  is  zero,  while  the  maximum  effect  on  rotation  occurs  not  only 
at  bottom  but  also  at  inflection  points.  Shear  force  or  bending  moment  also  experience 
changes  of  resonant  bands  in  the  same  fashion  as  that  of  the  rotation  or  displacement.  As 
expected,  the  maximum  effect  on  bending  moment  occurs  at  zero  slope  points  and  the 
maximum  effect  on  shear  force  occurs  at  inflection  points. 

In  addition  to  the  conclusions  made  in  this  study,  it  should  be  noted  that  the  effects 
by  coriolis  force  could  be  overestimated  because  downward  internal  flow  is  not  included 
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in  the  construction  of  a mathematical  model  even  though  the  downward  flow  tends  to 
reduce  coriolis  force. 

8.2  Further  Study 

In  this  dissertation,  Iwan-Blevin's  model,  which  is  one  of  the  self-excited  oscillator 
models,  has  been  implemented  as  a vortex  shedding  model  to  examine  the  effect  of  inter- 
nal flow.  In  that  model,  the  span  of  riser  is  divided  into  segments  that  are  within  the  reso- 
nance band  and  the  remainder,  which  is  outside  the  resonance  band.  However,  the 
numerical  application  of  the  spanwise  extent  of  resonance  band  is  not  feasible  due  to  the 
distortion  of  mode  shape  resulting  from  the  mixed  derivative  term  even  though  it  is  ex- 
pected to  be  small.  Thus,  we  need  to  develop  the  more  refined  current-vortex  model  to 
include  the  depth-dependent  resonance  described  above.  Besides,  the  nonlinear  coupling 
between  inline  and  transverse  vibration  and  the  extension  of  the  applicable  range  of  Re- 
ynolds number  may  be  required  for  the  development  of  the  refined  model  to  investigate 
the  effect  of  internal  flow  more  precisely. 

The  internal  flow  causes  the  distortion  of  trajectory  shape  as  discussed  in  the  pre- 
vious chapter.  Thus,  it  may  be  required  to  include  the  internal  torsion  in  the  riser  analy- 
sis. However,  the  coupling  problem  with  torsion  is  not  easy  task.  In  particular,  the 
torsional-axial  coupling  problems  still  remain  to  solve.  Uncoupled  internal  torsion  may 
be  applied  in  the  derivation  of  model  and  the  effect  of  internal  flow  on  it  may  be  investi- 
gated using  the  model  in  the  future  study. 

In  addition,  random  wave  excitation  and  corresponding  rig  motion  may  be  also  in- 
cluded in  the  future  study.  An  explicit  time-dependent  representation  of  wave  elevation 
and  velocity  may  be  generated  for  a time  domain. 


APPENDIX  A 

DERIVATION  OF  LINEAR  GOVERNING 
EQUATION  USING  HAMILTON'S  PRINCIPLE 


As  an  element  of  fluid  flows  with  a velocity  Vi  along  the  pipe  under  the  assump- 
tion of  small  deformation,  the  x -component  of  velocity  is  x + Vix7.  The  kinetic  energy 
in  a length  dz  of  the  pipe  is  that  of  the  pipe  plus  that  of  the  fluid,  which  is 

dT  = { \{mt  - mf)x2  + \mf[\2l  + (x  + Vix')2] } dz 

The  first  term  on  the  right  side  represents  the  kinetic  energy  of  the  pipe  and  added  mass, 
and  the  second  represents  the  kinetic  energy  of  the  internal  fluid.  The  strain  energy  due 
to  the  lateral  deflection  in  the  length  dz  is 

dV  =[\ EI(x")2  - \l(x')2]dz 


The  first  term  on  the  right  side  represents  the  stiffness  strain  energy  and  the  second  repre- 
sents the  tension  strain  energy.  The  negative  sign  in  front  of  tension  strain  energy  indi- 
cates that  it  reduces  the  total  strain  energy  due  to  the  lateral  deflection.  Hamilton's 
principle  states  that  5 {^(T  - V)dt  - 0. 

Substituting  the  kinetic  energy  and  the  strain  energy  and  performing  the  variation  gives 


{( mt  - m/)xdx  + ntf(x  + 


Vix/)(5x  + ViSx7)  - E\x" §x"  + Tx'Sx'  }dzdt 


Since 

3 

5*  = di 5x/  = 5 * " = 


each  term  may  be  integrated  by  parts  so  as  to  eliminate  the  various  derivatives  of  5x. 
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Performing  integration  by  parts,  we  obtain 

|2  f {mtx  + 2/w/Vi  x'  + mf\\x"  + (EIx")"  - (Tx')'  }6  xdzdt  = 0 
*t\  Jz 

where  all  the  integrated  terms  have  disappeared  because  of  the  boundary  conditions. 
Since  5x  is  arbitrary,  the  only  way  for  the  integral  to  be  zero  is  for  the  expression  in 
bracket  to  be  zero.  Equating  this  to  zero  gives  the  differential  equation. 

mtx  + 2/w/Vix7  + /w/V \x"  + (EIx")"  - (Tx')'  = 0 

Combining  the  external  hydrodynamic  force  and  the  lateral  component  of  the  difference 
of  the  hydrostatic  pressure,  we  obtain  Eq.  (2.68)  and  (2.69). 


APPENDIX  B 

DERIVATION  OF  FINITE  ELEMENT  MODEL  OF  LINEAR  SYSTEM 


To  produce  the  weak  form,  multiply  the  residual  error  function  R(x3)  by  a suffi- 
ciently smooth  test  function,  integrate  the  product  by  parts,  and  equate  the  result  to  zero. 
The  residual  error  function  of  Eq.  (3.54)  is  given  by 

R(x3)  = mix i + 2/w/Viii  + w/Vjxf  + (Elxf  )"  - wx^  - T tx'[ 

- a4pwD0uc(w  - xi)p(x3)  + \pwUcD0CDXi  (1  - p(x3)) 


and  test  function  v(x3)  = xi  (x3). 

The  variational  statement  is  given  by 

f Rx!fA:3  = [ [w,xi  + 2/w/ViXi  + mfV\x"  + (EIx")//  - vjx\-Tex"-a4pwD0ncpix^)w 

Jo  J0 

+ { jpwMc00CD(l  - P(x3))  + a4pwDouc  p(x3)}xi]xii&3  =0  (B.l) 


Integrating  terms  in  the  right  side  of  above  equation  by  parts  and  writing  them, 
j*  mjV\x'{  X\  etc?,  = mf\ \x\  xx  K - f mjV\x [ x[dx 3 

Jo  0 Jo 

f Tex"xi£A3  = TeX'jXil^  - [ Tex'x'^ 

Jo  u Jo 

(,(EIJt?)"xi<fc,  = (EIx")'Il|'n  - (EIx")5'  If,  + f (El x'px'ldx, 

Jo  u J° 


( B.2  ) 


Substituting  (B.2)  into  (B.l),  we  obtain  Eq.  (3.58),  i.e.,  the  weak  form  of  Eq.  (3.54).  Ap- 
plying same  procedure  for  Eq.  (3.55)  and  (3.56),  we  also  the  weak  forms  of  those  equa- 
tions as  seen  in  Eq.  (3.59)  and  (3.60). 
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Partitioning  whole  domain  into  a number  of  finite  elements,  Eq.  (B.l)  becomes 
n Cx3i+ 

E I [mtX\X\  + 2mfV\x\x\  - /w/V \x\x\  + ( E\x" )x'{  + (Tex^  )x\ 
m Jx3i 

+ {a4pwDoucp(x3)  + \pwucD0CD(\  - p(x3))  }xixi  ]dx3  + mfV\ x[x{  \ 1 

n f^3' 

+ {-(Tex{  )xx  +(EIx")/x1  -(EIx")x^) |*  = E I a4pwD0ucp(xi)wxidx3 

' 0 j=i  Jx 3/ 


Based  on  Galerkin's  approximation,  replace  above  problem  by  a finite  dimensional  prob- 
lem to  find  the  approximate  solutions  xih,  X2h,  and  Wh  e Hh such  that 

E '+[w,xihXih+ 2w/ViX/lhxlh-  W7/Vix/Ihxlh  + (EIx'/h)xi/h+ (Tex/lh)x/1 

w J*3/ 

+ {a4pwD0ucp(x3)  + \pw7jcD0CD(\  - p(x3))}xihXih]^3  + mf\ \x\hxx  | ^ 

+ {-(Tex'h )xlh  + ( El x" )' x lh - ( El x"  ) x'lh  |' 

P3' 

= E a4  pw  Do  wcp(x3)  vrh  x ih  <&3  ( B.3  ) 

w •’^Si- 


Approximated  solutions  and  test  functions  are  given  from  the  implementation  of  Hermite 
polynomial  as  basis  functions. 

4 

Xjh  (x3, 0 = E Xji  ( t)Ni  (x3  ) 

1=1 

4 

Xjh  (x3,  t)  = E Xji(t)Ni  (x3 

f=l 

4 

wh(x3,0  = E W|(0M(*3  ) j = 1,2  ( B.4  ) 

/=i 

The  substitution  of  Eq.  (B.4)  into  (B.3)  introduces  the  series  form  of  the  weak  form.  Fac- 
torizing the  coefficient  of  test  function  from  the  series  form,  we  have  the  compact  form 
such  as  Eq.  (3.44).  And  then,  we  can  obtain  the  matrix  dynamic  equilibrium  equation 
(3.61)  through  the  appropriate  choice  of  the  coefficients  of  test  function  as  shown  in 
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section  3.4.  Applying  same  procedures  for  Eq.  (3.59)  and  (3.60)  as  introduced  above,  we 
can  also  derive  the  corresponding  matrix  equation  (3.62)  and  (3.63). 
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